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These lecture notes discuss two effects which contribute to the reduction of the interference fringe 
contrast in matter interferometers. The first effect is the shot noise arising from a finite number of 
atoms used in experiments. Focusing on a single shot measurement we provide explicit calculations 
of the full distribution functions of the fringe contrast for the interference of either the coherent or 
the number states of atoms. Another mechanism of the suppression of the amplitude of interference 
fringes discussed in these lecture notes is the quantum and thermal fluctuations of the order pa- 
rameter in low dimensional condensates. We summarize recent theoretical and experimental studies 
demonstrating that suppression of the interference fringe contrast and its shot to shot variations can 
be used to study correlation functions within individual condensates. We also discuss full distribu- 
tion functions of the fringe amplitudes for one and two dimensional condensates and review their 
connection to high order correlation functions. We point out intriguing mathematical connections 
between the distribution functions of interference fringe amplitudes and several other problems in 
field theory, systems of correlated electrons, and statistical physics. 
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I. INTRODUCTION 
A. Interference experiments with cold atoms 

From the earliest days of quantum mechanics its probabilistic nature was the cause of many surprises and contro- 
versies pj. Perhaps the most unusual manifestation of the quantum uncertainty is a quantum noise: measurements 
performed on identical quantum mechanical systems can produce results which are different from one experimental run 
to another. At the level of a single particle quantum mechanics, the quantum noise is no longer a research topic but 
is discussed in undergraduate physics textbooks Q. However, the situation is different when we talk about quantum 
mechanics of many-body systems. One can ask seemingly simple questions to which there is no obvious answer: does 
it still make sense to talk about quantum noise when discussing measurements on many-body quantum states? How 
does the quantum noise manifest itself? Can one use this noise to extract nontrivial information about the system? 

The idea of the quantum noise analysis of many- par ticle systems is common to many areas of condensed matter 
physics [1, EL IE H| an d quantum optics @, fsL Ip. [Tol. fll L 



.12 L In the field of ultracold atoms it has been successfully 
16l . 17L [HI, |Ti| with many more theoretical proposals awaiting 



employed in a variety of recent experiments 

their turn (2(il . I2H I22I I23I HE HE Hy] . These lecture notes address a very specific experimental probe of the cold atomic 
ensembles - interference experiments. Our discussion focuses on a variety of interesting and important phenomena 
which originate from the fundamental quantum and/or thermal noise of cold atoms condensates and can be studied 
in interference experiments. Although focused on the specific type of experiments, the general methodology discussed 
in these lecture notes can be extended to a variety of other measurements on systems of cold atoms. 

Interference experiments constitute an important part of the modern toolbox for studying ultracold atoms. Original 
experiments used large three dimensional Bose-Einstein condensates (BEC) to demonstrate macroscopic coherence 
[27I . More recently interference experiments have been done with one and two dimensional condensates (28l . [2^ . l30l . 
l3ll. [32l. HH and demonstrated the important role of fluctuations in low dimensional systems. Matter interferometers 
using cold atoms [U, HE HE HE Hi, H3, HE S3, S3, E3, EE S3, S3, EE E3, EE El, have been considered 
for applications in accelerometry, gravitometry, search for quantum gravity, and many other areas (for a review see 
Ref. [51]). Interference experiments have been used to measure the condensate formation [EE Hj] as well as the 
critical properties of the BEC transition (54[. What is common to most interference experiments is the focus on 
the phase of interference patterns. Suppression of the fringe contrast is considered to be a spurious effect caused by 
noise and fluctuations. On the contrary, these lecture notes focus on understanding physical phenomena that underlie 
the imperfect visibility of interference fringes. As we discuss below, suppression of the fringe visibility comes from 
fundamental physical phenomena, such as the noise intrinsic to performing a classical measurement on a quantum 
mechanical wave function (shot noise) or classical and quantum fluctuations of the order parameter. 

In these lecture notes we discuss how one can use analysis of the contrast of interference fringes to learn about 
fluctuations of the order parameter [2E HE HE HE HE HI]- We will demonstrate that important information is 
contained not only in the average contrast but also in its shot to shot variations. For example, when we discuss 
fluctuating condensates, we will show that high moments of the interference fringe amplitudes contain information 
about high order correlation functions and thus provide valuable information about the system. 

The basic scheme of interference experiments is shown in Fig. [TJ Originally two condensates are located at distance 
d away from each other. At some point they are allowed to expand ballistically [59( until sizes of the clouds become 
much larger than the original separation between the clouds d. After the expansion the density is measured by 
shining a laser beam through the cloud. Interference leads to the appearance of the density modulation at a wave 
vector Q = md/Tit (see Fig. [T] and discussion in section HI)) . When the two condensates are coherent, the position of 
interference fringes is determined by the relative phase between the two clouds. Surprisingly the interference pattern 
will be observed even for two independent condensates which do not have a well defined relative phase (see e.g. Fig. [2]). 
To beginning readers it may seem confusing that we can observe interference in the absence of coherence between the 
two clouds. Or even more confusing, we discuss interference patterns when both clouds are number states and phases 
of individual condensates are not well defined. Several theoretical frameworks have been introduced to understand 
the origin of interference patterns in the absence of phase coherence [IE HE HE HE H3 • In this paper we explain 
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FIG. 1: Schematic view of the interference experiments with ultracold atoms. Two condensates, originally separated by the 
distance d, axe released from the traps and expand until the clouds overlap. The imaging beam measures the density of the 
atoms after the expansion. Quantum interference leads to the periodically modulated density projected on the screen. Projected 
density image is taken from the actual experimental data of Hadzibabic et al. 32] . 

1 mm 




FIG. 2: Interference pattern from the first set of experiments with sodium atoms. Taken from Ref. [27l |. 

the origin of interference fringes usingthe language of correlation functions and point out connections to Hanbury 
Brown and Twiss (HBT) experiments [7| in optics. This section provides a simple heuristic picture and a more formal 
discussion is left to section [Til 

What one measures in experiments is the density profile after the expansion p(r). Interference pattern appears 
as the density modulation p(r) = pQe~ ,( 3 r + c.c. + const. The absolute value of pq determines the amplitude of 
interference fringes and its phase defines the position of the fringe maxima and minima. Schematically we can write 
Pq ~ e l( " -i ^ 2 , where <j>\ t 2 are phases of the two condensates before the expansion (see discussion in section [TTJ) . In 
the absence of coherence (e 1 ^ 1-1 ^ 2 ) = 0, which implies that (/?q) = 0. Vanishing of the average, however, does 
not mean that interference fringes are absent in each individual shot. In the present case it only shows that the 
phase of interference fringes is random in each shot. We remind the readers that taking an expectation value in 
quantum mechanics implies averaging over many measurements. On the other hand we can focus on the amplitude 
of interference fringes and accept the fact that we can not predict their phase. Then we need to consider the quantity 
which does not vanish after averaging over the unpredictable phase difference. One such quantity is given by the 
density-density correlation function p(r)p(r') = |/OQ| 2 (e lC ^ r ~ r > + c.c.) + other terms. The right hand side of the last 
equation does not vanish when we average over random phases <^i,2 5 and we find finite expectation value 

(p(r)p(r')) = 2(| PQ | 2 ) cos(Q(r - r')) + const. (1) 

What this correlation function tells us is that in a single shot we can not predict whether at a given point r we will 
have a minimum or a maximum of the density modulation. However what we can say is that if there is a maximum 
at point r, it will be followed by another maximum a distance 2ir/Q away. While this simple argument explains the 



4 




FIG. 3: Schematic view of Hanbury Brown and Twiss noise correlation experiment as an example of intensity interferometry. 
Detectors at positions r, r measure the intensity of light coming from two distant incoherent sources. The "correlator" (denoted 
by the box) measures the coincidence events and thus the intensity-intensity correlation function. 

origin of interference patterns from independent condensates, it leaves many questions unanswered. For example, it 
is not obvious how accurately one can represent two independent condensates using states with a well defined but 
unknown phase difference. Also it is not clear how to generalize this analysis to elongated condensates, when we need 
to go beyond the single mode approximation and include phase fluctuations within individual systems. These lecture 
notes will present a uniform approach for addressing these and many other questions. 

When the focus of interference experiments is on measuring the phase, one usually averages interference patterns 
obtained in several shots. The result is easy to interpret: an average of many experimental runs is precisely what we 
define as a quantum mechanical average. However in experiments with independent condensates, summing interference 
patterns is not appropriate. The phase of interference patterns is random from shot to shot and adding individual 
images washes out interference fringes completely (for a nice experimental demonstration see Ref. (65j). Hence, in 
this case one needs to focus on interference patterns obtained in individual shots. In the absence of averaging, a 
single shot measurement contains noise. Thus to characterize such experiments we need to provide both the average 
value and the shot to shot fluctuations of the fringe contrast. The most comprehensive description of the fluctuating 
variable comes from providing its full distribution function. Theoretical calculations of the distribution functions of 
the fringe contrast will be the central part of these lecture notes. 

It is useful to point out the analogy between the approach discussed in this paper and the famous Hanbury Brown 
and Twiss experiments in optics. The original motivation for HBT experiments came from astronomy: the goal was 
to measure the angle between two incoherent stellar sources such as two different points on the surface of the star. 
Since the two sources are incoherent, this can not be done using a single detector: first order interference is absent 
and the measured signal is simply the sum of the two intensities [lH . The insight of HBT was to use two detectors 
and measure the correlation function of the two intensities as a function of the relative distance between the two 
detectors. One finds that this correlation function is given by 

(/(r)I(r')) ~ cos ((kr - k 2 )(r - r')) + const, (2) 

where ki 2 are wave vectors of photons arriving from the two points on the surface of the star, r, r' are positions of the 
detectors, and I(r), I(r') are the intensities measured in the two detectors (see Fig. [3]). Hence, the main idea of HBT 
experiments is that information is contained not only in the average signal (/(r)), but also in the noise. Such noise 
can be characterized by looking at higher order correlations. In astronomy, HBT experiments were used to measure 
several important properties of distant stars, including their angular sizes and the surface temperature [GfJ. 

There is an obvious analogy between equations (Q} and ([2]), but there is also an important difference. HBT 
stellar interferometers operate in real time, and averaging over time is built into the measuring procedure. In these 
experiments, fluctuations of /(r)J(r') — (I(r)/(r')), which would correspond to higher order correlations in I(r), are 
not easy to measure. On the other hand interference experiments with cold atoms are of a single-shot type: each 
measurement is destructive and gives a certain density profile p(j). A single image contains information not only 
about the two point correlation function (jo(r)p(r')), but about higher order correlations as well. The most important 
information for our purposes is contained in the interference pattern at wave vector Q. Essentially each interference 
pattern constitutes a classical measurement of the quantum mechanical operator pq . We remind the readers that we 
expect the phase of pq to be random, so the quantity of interest will be |/9q| 2 . By performing measurements several 
times we will find not only the average value of this operator, but also its higher moments. Ultimately we should 
be able to reconstruct the entire distribution function for |/9q| 2 . So the simplified and idealized procedure that we 
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analyze is the following: one performs interference experiments many times. Each experiment is analyzed by doing 
a Fourier transform of the density to extract pq. The histogram of the measured values of \pq\ 2 will be the main 
subject of these lecture notes. We will demonstrate the wealth of information that can be extracted from analysis of 
such histograms. 

As a passing note, we mention that the setup considered in Fig. [1] is not the only possible configuration for 
interference experiments. In another common setu p o ne makes several copies of the same cloud using Bragg pulses, 
and observes an interference between them [63, [68ll69l l70l . l7ll . l72l . l73j . 

It is useful to put our work in the general perspective of noise analysis in physical systems. Understanding photon 
fluctuations is at the heart of modern quantum optics and provides a basis for creation, detection, and manipulation of 
non-classical states of light [ll[ . The field of quantum optics has a long and fruitful tradition of using the higher-order 
correlation functions as well as the shot noise to characterize the quantum states of light. The notion of higher-order 
degree of coherence was first introduced by R. Glauber in 1963 Q, also by Klauder and Sudarshan Q and by Mandel 
and Wolf [l(| ■ The knowledge of photoelectron counting distribution function reveals such non-classical features 
of light as antibunching (jij, sub-Poissonian statistics and probe of violation of Bell inequalities. In particular, 
third-order correlations provide a test for distinguishing between quantum and hidden variable theories in a way 
analogous to that provided by the Greenberger-Horne-Zeilinger test of local hidden variable theories (75j . Interference 
of independent laser beams was first observed in Ref. [76J and stimulated a number of theoretical studies (for reviews 
seeRefs. [l3, IzE Iz3 ) ■ 

In condensed matter physics, noise analysis was also suggested as a powerful approach for analyzing electron systems 
Q. It was demonstrated theoretically that in certain mesoscopic systems current fluctuations should contain more 
information than the average current itself. In particular, the third and higher moments contain important quantum 
information on interaction effects, entanglement and relaxation processes (see e.g. Refs. fZ^HO])- Specific proposals 
exist for detecting statistics of quasi-particles [8l| , understanding transmission properties of small conductors [j] and 
observing entanglement between electrons Q • Perhaps the most spectacular experimental success of the noise analysis 
in electronic systems has been the demonstration of the fractional charge of quasiparticles in the fractional quantum 
Hall regime [6j . However generally the noise analysis in condensed matter systems did not become the detection tool 
of the same prominence as in quantum optics. The main reason for this is the excruciating difficulty of the noise 
measurements in solid state experiments. One often needs to measure a signal which is only a part in a million of the 
unwanted technical noise. 

In the field of ultracold atoms experiments analyzing quantum noise are only starting. However we have already seen 
spectacular successes in several recent experiments. Analysis of noise correlations in the time of flight experiments 



20] was used to demonstrate fermionic pairing [13[ as well as HBT type correlations for atoms in optical lattices 
14 ITU [l6j](see also I. Bloch's lectures in this school). Single atom detectors have been used to demonstrate HBT 



noise correlations for cold atoms [Tt], EH E^l • Strongly interacting systems of cold atoms are expected to realize 
analogues of important models of condensed matter systems [12, HE H3| • Being able to study noise in such systems 
should provide an important new perspective on strongly correlated states of matter and have a profound effect on 
many areas of physics. We hope that these lecture notes will stimulate more experimental work in analyzing noise in 
interference experiments. The first success in this direction was the recent observation of the Berezinskii-Kosterlitz- 
Thouless (BKT) transition [H,[86| by Hadzibabic et al. in Ref. [II]. 



B. Fundamental sources of noise in interference experiments with matter 

Two fundamental sources of fluctuations in the amplitude of interference fringes are the shot noise and the order 
parameter fluctuations within individual condensates. 

Shot noise comes from the finite number of atoms used in the experiments. Let us discuss limiting cases first. 
Consider an interference experiment with one atom. Before the expansion the atom is in a perfect superposition of 
being between the two wells. After the expansion we get a perfectly periodic wave function ip(r) = 2C cos(^±^) 
(for a more detailed discussion see section ITT1) . The expectation value of the density operator is (p(r)) = (|V>( r )| 2 ) = 
2|C| 2 (cos(Qr + </>) + 1). However this average value will not be measured in a single shot. A single measurement finds 
the atom at a single point. The expectation value of the density determines probabilities with which we can find atom 
at any given point, but in a single measurement we collapse a quantum mechanical wave function and observe the atom 
at one point only. Can one reconstruct the entire amplitude of the interference pattern pq = \C\ 2 e~ z '> > from a single 
measurement? Obviously the answer is no. In the opposite case of a very large number of atoms in the same single 
particle state one should be able to reconstruct a complete interference pattern already from a single measurement, 
since the measurement of positions of many atoms performs a statistical averaging implicit in quantum mechanics. 
In the general case of experiments with a finite number of atoms, the question arises how well one can determine the 
amplitude of interference fringes from doing a single shot measurement. Formulated more accurately the problem is 
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FIG. 4: Simplified setup of interference experiments with ID Bose liquids (see e.g. Ref. [29|]). Two parallel condensates are 
extended in the x direction. After atoms are released from the trap, clouds are imaged by the laser beam propagating along 
the z axis. Meandering structure of the interference pattern arises from phase fluctuations along the condensates. The net 
interference amplitude Pq(L) is defined from the density integrated along the section of length L. 



to determine probabilities of finding a certain amplitude of interference fringes, |pq| 2 , in a single measurement. 

Fluctuations of the order parameter are particularly important for low dimensional systems. If the condensates are 
confined in one [83, [8^] or two [3l|, [32, [H, lS3, Ha] dimensions, then the true long-range order may not exist. Rigorous 
theorems forbid true lon g ra nge order in two dimensional systems at finite temperature and in one dimension even 
at zero temperature (90l l9ll. |92 |. What this means is that low dimensional condensates can not be characterized 
by a single phase and we need to take into account spatial fluctuations of the order parameter. Effects of such 
fluctuations on interference experiments are illustrated schematically in Fig. 01 Two one dimensional clouds expand 
in the transverse direction. Each point along the condensates has a local interference pattern, but in the presence 
of phase fluctuations (either thermal or quantum), these patterns are not in phase with each other. It is natural to 
define the net interference amplitude from the density integrated over the axis of the system (the so called columnar 
density). In many experiments such integration is done by the measurement procedure itself. For example, systems 
such as shown in Fig. [4] originally had imaging done along the axis of the interferometer [29]. Then the laser beam 
integrates the atomic densities within the imaging length. Integrating over local interference patterns which are not 
in phase with each other leads to a reduced contrast of the net interference fringes. 

In earlier literature smearing of interference patterns by fluctuations was considered an unwanted effect (27j . The 
point of view presented in these notes is quite the opposite. Suppression of the fringe contrast is an interesting effect 
which tells us about important phenomena in Bose condensates, such as thermal and quantum fluctuations of the order 
parameter. By analyzing such suppression we can extract non-trivial information about the system. In particular, 
it has been shown in Ref. [55[ that the scaling of the average interference signal with the observation area contains 
information about the two-point correlation functions within each cloud. Recent experiments [32j by Hadzibabic et 
al. used this approach to observe the BKT transition in two dimensional condensates (see discussion in section IIVI) . 
We also note that such experiments can be used to extract information which is difficult to obtain by other means. 
For example, in section IIVI we discuss that analysis of the high moments of the contrast tells us about high order 
correlations within individual clouds. 

These lecture notes are organized as follows. In section [TT] we discuss how interference fringes appear for ideal non 
interacting 3D BEC's at zero temperature. In section ITO1 we analyze the shot noise for ideal condensates. In these 
lecture notes ideal condensates are understood as clouds of non-interacting atoms which before the expansion occupy 
a single mode within each of the traps. The problem of interference of independent condensates of ideal bosons has 
been extensively analyzed in the literature before [13, [fill) Is2) IsS [H, [93|, 0, [9f| [9(| [93, HI]. In particular in an 
important recent paper [99| , Polkovnikov showed that the variance of the fringe amplitude decreases as an inverse 
power of N, with a non-universal coefficient which contains information about the state of each cloud (e.g. coherent 
states vs Fock states). In this paper we develop a general formalism for calculating the full distribution functions of 
the fringe amplitudes in interference experiments with ideal condensates with a finite number of atoms. We apply this 
formalism to obtain distribution functions for several experimentally relevant cases such as states with a well defined 
phase difference between the two clouds and Fock states of atoms. Effects of the order parameter fluctuations are 
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discussed in section HVl We obtain distribution functions for both one and two dimensional condensates in the limit 
when the number of particles is large and the shot noise can be neglected. We also discuss intriguing mathematical 
connections between these distribution functions and a number of impo rtant problems in physics, such as the quantum 
impurity p robl em in a low dimensional interacting electron system [lOOl j or the distribution of roughness in systems with 
1/ f noise |10l| . In these lecture notes we do not address the issue of technical noise which is obviously important for 
understanding real experiments. In the concluding section|Vl however, we comment on the experimental requirements 
for observing some of the phenomena discussed in these lecture notes. 



II. INTERFERENCE OF IDEAL CONDENSATES 



In this section we discuss why interference fringes appear for ideal non interacting BEC's at zero temperature. We 
follow Ref. [Hj], and introduce notations for subsequent sections. First we consider the case of two clouds with a well 
defined relative phase, where appearance of interference fringes can be understood at a single particle level. Then we 
show that almost ideal interference fringes appear even when two expanding clouds are uncorrelated, provided that 
the number of particles in each cloud is large. 



A. Interference of condensates with a well denned relative phase 



1. Basics of interference experiments. First quantized representation 



To illustrate how interference fringes arise, let us start by considering a simple case of two BEC clouds with a well 
defined relative phase. Here we neglect interactions between atoms, so initially all atoms are assumed to be in the 
same single-particle state (single mode approximation). After the confining potential is removed, the single-particle 
state evolves with time, but many body wave function remains in the product state. The interference appears as a 
result of single particle wave function evolution, which can be studied in detail. 

Normalized single particle wave functions for two clouds will be denoted as ipi(r, t) and ip2(^,t), and the initial 
relative phase is (p. If the total number of particles equals N, then the complete wave function of the system in the 
first quantized notations at any moment of time is given by 

N 1 

*(r 1; rjv , t) = J] -=0Mr„, t)e^' 2 + </> 2 (r n , t)e-^' 2 ). (3) 

n— 1 * 

This wave function satisfies the proper symmetry requirements for permutations of and Tj , and evolution of ipi (r, t) 
and ^(r, i) is controlled by the single particle Schrodinger equation. Initial overlap of the states ipi( r ) = V'i( r ; 0) and 
■02 (r) = ^(r, 0) is assumed to be negligible: 

'dr^J(r)^a(r)wO. (4) 

The expectation value of the total density corresponding to the wave function © is 

(p(r,i)) = ^ (|^i(r,t)| 2 + |^ 2 (r,t)| 2 + 2Re [ e ^i(r,i)^(r,t)]) . (5) 

The expectation value of the total density displays an interference pattern due to the last term of Eq. ©. As a 
simple example, let us assume that tpi( r > *) an d ip2(r,t) are initially in the Gaussian states centered at points ±d/2, 
and their widths are Rq <C d. Then the evolution of the single particle wave functions can be simply calculated, and 
the result is 

^ (r-d/2) 2 (l + ifit/ m fig) 

, _ (r + d/2) 2 (l + ifit/mH2 ) 

where the widths of the wave packets, Rt, at time t are given by 
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FIG. 5: Schematic view of the interference experiment with 3D condensates. 

We will be interested in the regime, when the sizes of the clouds Rt are much larger than original distance between 
the clouds, that is 

i? t >d>i? . (9) 

In this regime, the clouds overlap strongly, and the real parts of the exponents in Eqs. (O and ([7]) are responsible for 
the broad overall density profile. Imaginary parts in the same exponents give rise to interference effects in the last 
term of Eq. ((5]). Thus the interference part of the density is equal to 

N - r2 +f /4 fh rd \ 

MlF 1 ' C0S UiPF^| (10) 

For sufficiently large t, one can substitute Rt ~ ht/mRo, and obtain oscillations of the density at wave vector 
Q = md/fvt, with positions of the minima and the maxima controlled by the relative phase (p. The Fourier transform 
of the density at wave vector Q is 

Physically Q can be understood as the momentum difference of the two particles which have been released from 
the two traps and arrive simultaneously at the detection point r. This can be seen from the following quasiclassical 
argument. A particle released from the condensate one and detected at time t at point r has momentum 

Qi =m(r-d/2)/(fii). (12) 

During the expansion this particle picks up a phase Qi(r — d/2). Analogously a particle originating from the conden- 
sate two has momentum 

Q 2 = m(r + d/2)/(fit) (13) 

and picks up a phase Q2(r + d/2). The interference pattern arises from the oscillating structure of the phase difference 
with the wave vector of oscillations 

Q = Qi Q2 = md/ht. (14) 

This simple argument shows that as long as the original sizes of the clouds are much smaller than the distance between 
them, after sufficient expansion one should observe oscillations of the density at the wave vector Q determined by the 
distance between the clouds. 

The interference patterns which we introduced up to this point appear as a result of the time evolution of single 
particle states. The many-body nature of the state comes into Eq. ([5]) only as a prefactor N. In principle, one could 
have done the same experiment with only one particle. The same result as Eq. ([5]) can be obtained in this case by 
doing experiments many times and averaging over individual experiments: in each particular realization the particle 
is observed at some random point r with probability (p(r,t)). In experiments with a large number of atoms, N, used 
in each shot, each absorption image is a result of N measurements of single particle wave functions. This performs 
the statistical averaging implicit in quantum mechanics, and leads to the density profile close to Eq. ^ in each shot. 
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2. Second quantized representation 



To set up the stage for later we will now present the discussion leading to Eq. (fTTjl using the second-quantized 
formalism. Wave function §3§ in the second quantization at t — can be written as 



^N) = w ^(^ /2 + -h^ /2 ) N \o). 



Here a\ and a 2 are creation operators for clouds one and two: 



dripi^tp' (r) 



ip(r) is the second quantized operator for the boson field, which satisfies the usual commutation relations 



i>(r'),^(r) = <5(r'-r), ^(r'),V>(r) = ft(r'),ft(r) 



= 0. 



Operators a\, a 2 and their conjugates satisfy the canonical boson commutation relations: 



Sij, [<H,a,j] 



at, a] 



= 0. 



(15) 



(16) 



(17) 



(18) 



Different initial states can be simply written using the Fock basis of operators a\ and a\. For example, the initial 
state for two independent condensates with N\ and N 2 particles in clouds 1 and 2 can be conveniently written as 



\N U N 2 



1 



VNx\N 2 \ 



(ai)*v 2 no>. 



(19) 



Initial state written in the Fock basis of operators a\ and a 2 contains all information about properties of the interference 
amplitudes. During the free expansion the occupation numbers of states one and two do not change, and only the 
single particle wave functions ip\{r,t) and ip 2 (r,t) evolve. After the expansion, the many-body wave function at time 
t can be obtained from the initial state written in the Fock basis of operators a\ and a 2 using substitutions 



a\ — > / drijj 2 (r,t)$ '( r )- 



For example, substituting Eqs. (|2U1) - (|2T|) into Eq. (fTS")) . the wave function ([3]) considered earlier is written as 



\ip,N,t) 



dr(^i(r,i)e iV / 2 + 7/; 2 (r,t) e -^/ 2 )^(r)t |0) 



N 



(2 N N\y/ 2 

In the long time limit ([9]) considered earlier, single-particle wave functions ipi(r, t), ^2(1", t) can be written as 

V>iM) = Ul (MK Qir , 



(20) 
(21) 

(22) 



(23) 
(24) 



where Qi, Q2 are defined by Eqs. (|Tl2|) - (fT"3"f and u\{r, t), u 2 (r, t) are slowly varying real functions, which determine the 
overall density profiles. Since clouds overlap strongly after the expansion, these functions are normalized according 
to 



J wi(r,i) 2 dr = 1, / u 2 (r,t) 2 dr = 1, 
J u 1 (r,t)u 2 (r,t)dr w 1. 



(25) 
(26) 



The operator, which corresponds to the amplitude of density oscillation at wave vector Q is written in the second 
quantized notations as 



p Q = J dvp(v)e^ = J , 



iQr 



(27) 



10 



To find out the statistical average of the amplitude of density oscillations for state (|22[) . we need to evaluate the 
following matrix element: 

(pq) = (<p,N,t\fa\<p,N,t) = (<p,N,t\ I dr^(r)^(r)e^\^ 7 N 7 t). (28) 



To evaluate such matrix elements, first we need to know how annihilation operator ip(r) acts on a state \ip, N, t). Since 
\ip,N,t) is obtained from \ip,N) by substitutions (|20 p -(|2i p . it is easy to see that 

{tp,N,t\ft(T)$(r)\tp,N,t) = {(p,N\ (ot« 1 (r,*)e- << ' ir + 4«2(r > *)e~ <q ") (a lUl (r,t)e l ^ r + a 2 u 2 (r,t)e^ r ) \ip,N), 

where \<p,N) is defined in Eq. (| 1 5|) and is written only in terms of a\ and a\. Integration over dr in (|28[) can be 
done using normalization (|26p . and assuming that Ui(r,t) vary at scales much larger than l/Q. Since Q = Qi — Q2, 
evaluation of Fourier transform picks only the product of the first term in the first parentheses and of the second term 
in the second parentheses of the equation above. We obtain 

(Pq) = (<f,N\a\a 2 J ui(r, t)u 2 (r, t)dr\ip,N) = 

&,N\a\a 2 \p,N) = ( V ,N-l\(J^e- i v/* ) J \<p, N —1) = ^-e~ iip , (29) 

which is the same as obtained from a single particle discussion in Eq. (|11| . 

The example above illustrates how the matrix elements of many-particle operators at time t can be evaluated using 
initial states written in the Fock basis of operators a\ and a\. In general, when one needs to evaluate an expectation 
value of some normal ordered combination of operators ip(r), and xf. >>(t) over the final state at time t, one needs to 
make substitutions 

i>(r) -> aiui(r, t)e lQir + a 2 u 2 (r, t)e^ r , (30) 
V3 f (r) -» a\ Ul {r, t)e~ zQir + a\u 2 (r, i)e" lQ2r , (31) 

and evaluate matrix elements over the t = state, written in the Fock basis of operators a\ and a 2 . It is important that 
the expression needs to be normal ordered using commutation relations (|17|) before making substitutions ()30|) - ()31l) . 
since after substitutions (j30)) - ([3"Tj) fields ip(r), ijj >(r) do not satisfy the exact commutation relations (|17p . 

Another way of seeing this is to realize that in Heisenberg representation (see e.g. chapter 6 of Ref. |l02j |) 
substitutions (|30p -()3i p perform the time evolution of a product of boson annihilation operators 

^{r x ,t).4{v n ,t) = e lflt ^{v 1 )...^{r n )e- lflt (32) 
only when this product acts on states of the form 

|*o) = ^^(4)^(4)^10). (33) 

Thus for calculating the expectation values of the form 

{^\e i6t p{v 1 )...p{v n )e- i6t \^) (34) 
one first needs to normal order p(r ■{)... p(r n ) using commutation relations (|17p as 

p(ri),..p{r n ) = /m(ri, ...,r„)-0 t (ri)...-0 t (r m )V'(r m )...V'(ri), (35) 

m < n 

and only then use substitutions (p?U|) -(p?T j) to evaluate matrix elements: 

(* \e i6t p(r 1 )...p(v n )e- i " t \* ) = £ / m (n, r n )<* |e iA V t (ri)»^ t (rr B )e- iA V^(r m )..^(r 1 ) e - <At |*o> = 

711 < n 

^/ ro ( ri ,..,r n )(^ | (almin,^^ +4 U2 (r 1) t)e- i ^ ri ) ... (a\ Ul (r m , t) e -^™ + 4u 2 (r m , t) e -^^) x 

(a lUl (r m ,t)e iQ ^ + a 2 u 2 (r m ,t)e iCl ^) ... (a lUl (n, t)e^ + 031*2 (ri, i)e~ jQ2ri ) |* )<36) 
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B. Interference of independent clouds 



The surprising phenomenon which was observed in Ref. [27| is the appearance of interference fringes in the case 
when condensates are completely independent. To illustrate how interference fringes appear for independent clouds, 
let us now discuss the case when the numbers of particles in each of the clouds, Ni and N 2 , are fixed, hence the phase 
difference between the two clouds is not well defined. Initial state in the Fock basis in this case is given by (fT9|) : 



l*x,* a > = -^(«t)*(4)*l<)>. (37) 



Using the formalism of second quantization explained earlier, one can evaluate (pq) by analogy to Eq. (|29|) : 

(pq) = (N u N 2 \a\a 2 \N 1 ,N 2 ) = 0. (38) 

However, (pq) = doesn't imply that there are no interference effects for independent condensates. Indeed, (pq) 
gives only the statistical average over many experiments, according to the usual interpretation of expectation values 
of operators in quantum mechanics. Being a quantum operator, pq has non vanishing quantum fluctuations. In each 
particular realization of experiment, complex number pq can have a nonzero value. To show this, let us evaluate 
(I/°q| 2 )' which is the density-density correlation function at wave vector Q : 



\PQ\ 2 ) = (PQP-Q) = ( / drdr'V-t( r )^( r )v;t( r ')^ (r ') e iQ(r-r') ) = 

( / rfrdr'-0t( r )^t( r ')^( r )^( r ') e <Q(r- r ') + f dr^(r)i>(r)). (39) 



These matrix elements can be evaluated using the second quantization prescription of the previous section, and the 
result is 

(\pq\ 2 ) = (N u N 2 \a\a\axa 2 + a\ax + a\a 2 \Nx,N 2 ) = N X N 2 + AT X + N 2 . (40) 

In the limit of large N\ = N 2 — N/2, the leading contribution to (|pq| 2 ) is the same as for the state with the fixed 
phase. Information about the full distribution of the quantum operator pq is contained in higher moments of the 
distribution. If one considers higher moments of the type (\pq\ 2ti ) = (PqP-q): the leading contribution in the limit 
of large Ni and N 2 will again have the form 

(Ipq| 2 "> = (JW (l + 0(i-) + 0(i-)) . (41) 

Corrections which appear because of the normal ordering result in sublcading terms which are denoted by 0(1/ Ni) + 
0(1/N 2 ). The leading term implies that in the limit of large N\ and N 2 , the distribution function of |pq| 2 is highly 
peaked near the value N±N 2 , with the relative width which is proportional to the inverse square root of number of 
particles. Any operator of the form PqP™q will have zero expectation value for m ^ n similar to pq, which means 
that the phase of the complex number pq is uniformly distributed from to 2ir. The expectation value of any operator 
which depends on the phase of pq becomes zero due to the averaging over the phase. 

The physical picture which emerges from the calculations is the following [6l], [13, EH : for two independent ideal 
clouds in the limit of large N the absolute value of interference fringe amplitude is the same as for the state with 
a fixed relative phase, but the position of the intensity minima fluctuates from shot to shot. The state with a fixed 
number of particles is a superposition of states with fixed relative phases. For example, 

W2,iY/2)=(f) 1/Z XS'^' (42) 

In the limit of large N the phase states are almost orthogonal, and the measurement picks some value of the relative 
phase. Since the relative phase is not well defined for independent clouds, in each particular experiment the positions 
of the minima will fluctuate from shot to shot. To distinguish independent clouds from states which have correlated 
relative phases, one needs to do a series of experiments and measure not only the absolute magnitude of interference 
fringes, but also the positions of the minima. Experiments which distinguish states with a fixed relative phase from 
some other many body states are already being done, and can be used i.e. to measure the temperature [lOSj or to 
study the dynamical evolution [2^, [3(1 EH, EE Ea| of the relative phase. 
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III. FULL COUNTING STATISTICS OF SHOT NOISE 



As has been explained in the previous section, for experiments with independent clouds the average interference 
amplitude depends only on the number of particles per cloud. In this section we consider not only the average 
interference amplitude, but also its shot to shot fluctuations due to a finite number of atoms in the clouds. We 
will demonstrate that while the average value of |pq| 2 depends only on the number of particles per cloud, the full 
distribution function of the variable 

R=\ PQ \ 2 (43) 

contains information about the states of individual clouds. Our analysis is motivated by the earlier work of Polkovnikov 
[9^ | , who showed that the variance of the fringe amplitude decreases as the inverse power of the number of particles 
per cloud, N, with a non-universal coefficient which contains information about the state of the clouds. Experimental 
observation of the effects discussed in this section requires systems with a small number of atoms. This may be 
realized with micro-BEC's on chips 0, M, EI ES E3, El E! tHO] • 

The shot noise for finite N has a fundamental nature, which stems from the probabilistic nature of quantum 
mechanics. Distributions of R obtained below correspond to the following "idealized" experimental procedure: release 
the confining potential and take an absorption image of the columnar density on an ideal CCD camera with 100% 
efficiency (photon shot noise is ignored). To obtain the amplitude of interference fringes, pq, extract a Fourier 
component of the density at wave vector Q from each image separately. The results of many experiments give the 
histogram W(R) of the values of R = |pq| 2 - We note that the quantum observable R defined in such way is a many- 
body operator, and calculation of its full distribution is a non-trivial task, even when all atoms are in the same state, 
such as for the case of a well defined relative phase between atoms in the two wells. In this section we develop a 
general method to find distribution functions of R analytically. 

We note that in our idealized setting we find interference patterns at a well defined wave vector Q. We expect that 
the finite size of the systems in transverse direction after expansion and collisions during the initial stage of expansion 



broaden the peak in the Fourier space to a finite, but small range of wave vectors around Q 104]. Hence a one should 
consider R = \pq\ 2 as an integral over the peak in the Fourier image of the density. Also in this paper we will discuss 
the amplitude of the interference fringes whereas experimental papers typically discuss visibility of the interference 
patterns. The two quantities differ only by the trivial rescaling. 

Results for independent clouds in coherent (solid) and number (dashed) states for N\ = N2 = 100 are presented 
in Fig. [Sj One can see, that even for a relatively large number of atoms, N = 100, fluctuations due to shot noise 
are appreciable . In Fig. [7] we compare the full counting statistics of R — \pq\ 2 for independent clouds in coherent 
(dashed) and number (dotted) states with N\ = N2 = 20, and for clouds with a fixed relative phase (solid) with total 
number of atoms N = N\ + N2 — 40. Distribution functions for the cases of (i) well defined relative phase between the 
clouds and (ii) fixed number of atoms in each cloud are very close. They become indistinguishable when R is rescaled 
by its average value (i?), although each of them differs considerably from the Gaussian distribution. The distribution 
function of R is wider for coherent states compared to number states, as was suggested in Ref. [99| based on the 
study of the variance of the two distributions. Hence the conclusion is that coherent and number states can be easily 
distinguished based on the statistics of fluctuations of R relative to its average value. 

In principle, when the two clouds are prepared with the same relative phase over many experiments, it is possible 
to distinguish independent condensates in number states from states with a fixed relative phase using a set of several 
interference experiments: minima positions arc uniformly distributed for independent clouds, while for states with 
a fixed relative phase positions of the interference minima are always at the same points in space. However one 
can imagine the situation when the clouds are prepared in state with a fixed relative phase but the relative phase 
itself is random from realization to realization. Our results show that it is practically impossible to distinguish such 
states from number states by looking at the distribution of the amplitude of interference fringes. Our calculations 
below provide additional support to the physical interpretation [oil . l62l . [63], [64| presented in the previous section. 
Fluctuations of the absolute value of the interference amplitude are the same for the cases when clouds have random 
relative phase and when clouds are prepared in number states so the random relative phase is " measurement induced" . 

The method developed in this section can be generalized to a variety of experimental situations, i.e. several 
independent condensates. Different squeezed states within individual condensates can be considered, and measurement 
of full counting statistics of shot noise can be used as an experimental probe to distinguish between different correlated 
states. 

As noted earlier, we will be interested in the full distribution function of the positive definite quantum observable 

R = \p Q \ 2 = PQP-Q, (44) 

defined by Eq. (|3"5|) . To calculate its full distribution function, W(R), one needs to know expressions for higher 
moments (R n ). After that, one has to solve the "problem of moments", i.e. to recover the distribution function on 
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FIG. 6: Rescaled distribution functions of R — |pq[ 2 for independent clouds in the coherent states (solid) or in the states with 
a well defined numbers of atoms (dashed). Here TVi = N2 = 100. 
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FIG. 7: Distribution functions of R = |pq| 2 for independent clouds in the coherent (dashed) and in the number (dotted) 
states with Ni = N2 — 20. Solid line is a distribution function of R for clouds with a fixed relative phase with total number 
N — Ni + N2 = 40. Distribution functions for states with a fixed relative phase and with fixed numbers are very close, and 
become indistinguishable when R is normalized by its average value {R). 



the (0, 00) interval using all moments. In general, this procedure is numerically hard and unstable, unless higher 
moments have a certain analytical form. If the expression for higher moments is known analytically, then one can 
sometimes avoid the "problem of moments" by calculating the so-called characteristic function, xWi which is the 
Laplace transform of W(R) : 



X(A)=/ e-^W{R)dR= E ±^P-W(R)dR = £ ( ~ A) ^ } (45) 



i=0 i=0 

If x(A) can be calculated analytically, then W(R) can be recovered by the inverse Laplace transform. In o ur case it 
is more practical to calculate not the characteristic function, but the analog of the Hankel transformation jl05l . Il06j | 
of W(R), given by 

- E w ( ^ l> - (46) 

n=Q V 

Using the expansion of the zeroth order Bessel function, one can write 



Z(iX) = I W(R)J Q (2XVR)dR. (47) 

(48) 



The inversion of the transformation can be found using the orthogonality condition for the zeroth order Bessel functions 
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/ °° J (Aa;)Jo(Ay)|x|AdA = S(\x\ - \y\), which gi\ 

W(R) = 2 / Z{iX)J (2XVR)XdX. (49) 



By the end of this section we will provide analytical expressions for Z{iX) for certain cases (see Eqs. (|62|) . (|69)l . ([T2])). 
from which T4 7 (i?) can be obtained by simple numerical integration according to Eq. (|4"9"1) . 
To proceed we note that /5q and P-q commute with each other: 

[PQ,P-q]=0. (50) 

Operators /5q are understood as in Eq. ([27)) . without the projection on single particle states tpi(r,t),ip2(^,t) as in 
Eqs. ([30| - ([3T|) . The latter substitutions can be only done after normal ordering (see discussion in section Hi A 2 [) . 
Hence we find 

Z (i\) = ( r ^Kp^+p-^)). (51) 
Jo 2?r 

Indeed, after expanding the exponent and integration, only even degrees of iX survive, and non vanishing terms are 
exactly what is needed for Eq. (|46|) . 

The normal ordering of Eq. (|51|) can be done using the following identity: 

e //(r)Vit(r)^( r )dr = . e JV< r > - 1)</>+ (r)^(r)dr . _ ^) 

Here we have assumed that operators ?/^(r), ^ (r), have the canonical commutation relations given by Eq.JT7|). Eq. 
(|52| is a generalization of the simpler identity [l07| | : 

e Aata = . e (e A -l)at a . ^ 

for operators which obey the commutation relations 

[0,0^ = 1. (54) 

The normal ordering signs : : mean that all creation operators should be put to the left of annihilation operators in 
Taylor expansion of expressions being ordered. To illustrate the meaning of Eq. (|53|) . let us consider the expansions 
of left and right side up to A 2 . The left hand side is 

\2 \2 \2 

1 + Aa f a + yaWa + 0(A 3 ) = 1 + (A + y )a f a + yaVaa + 0(A 3 ). (55) 
The right hand side is 

,A 1 \2 



1 + (e A - 1) : a) a : +^ ^- : a) aa) a : +0(X 3 ) 



2 

1 + (A + y + 0(A 3 )) : at a : + ^!±^!l . a t aa t a . +0(A 3) = 

A 2 A 2 
1 + (A + y ) flta + Y atfltaa + 0(A 3 ). 

Eq. (|53[) holds not only up to A 2 , but to all orders in A and plays an important role in quantum optics. 

Using the definition of pq given by Eq. ([2T)l . one can apply Eq. (f5"2"|) with /(r) = 2iA cos (Qr + 1^9) and rewrite 
Z(iX) in Eq. dH]) as 

Z (iX) = (: ^ e /{e ,U » (t! ' + "-^ , (r)te<r (56) 

Jo 27r 
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A. Interference of two independent coherent condensates 



Let us first explain how to evaluate Z(i\) for independent clouds in coherent states [8J, 11081 ] of operators a\ and a 2 
with eigenvalues V^/V xe 1 ^ 1 and ^/N 2 e 1 ^ 2 ■ Since coherent states form a complete basis, any initial state can be expanded 
in this basis, and thus the problem of calculating of W(R) is essentially solved for arbitrary initial states. Coherent 
states are convenient, since they are the eigen states of the annihilation operator, and the annihilation operator acts 
on them as a c— number. Hence after making substitutions (|30|) - <|3 1 [) into the normal ordered expression, one can 
substitute operators a,,aj by numbers V^/Vje 1 ^*, \fN i e~ l ^ i . Since the normal ordered expression is obtained by the 
normal ordering of the Taylor expansion of Eq. (|56|) . one needs to collect the Taylor series back. For coherent states, 
the whole procedure is equivalent to removing the normal ordering signs in Eq. (|56p . making substitutions (|30|) - (|3 1[) 
and treating operators Oj, a\ as numbers V^V^e*^*, VN i e~ % ^ i . Thus we obtain 

Z(iX; y/Nie^ 1 , VN 2 e^ 2 ) = ^ e /( e2iAcos(Qr+ ^-0(^i"i( r -*) 2 +^ M 2( r >^ 
1 ' Jo 2?r 

Similar to section Til A 21 we assume the that normalized functions u±(r,t) and u 2 (r,t) strongly overlap and vary at 
scales much larger than 1/Q, which is equivalent to 

J e mQr u a {r,t)up(r,t)dr = S n0 . (58) 

Then integration over dr in the exponent of Eq. (|57p can be done using the following equations: 

^ 2l Acos(Qr+ v ) _ A (7VlWl (r,t) 2 + N 2Ul (r,t) 2 )dr = 



"■~° c (A \\2m f 

E -7^T /(^Vi«i(r,t) 2 +iV 2 u 1 (r,t) 2 )(e^ r +^+ e - l (Q r +^) 2m dr = 



(Ni + N 2 ) [t > = ( Jq(2A) _ 1)(Ni + Na) m 

"—^ {2m)\ rrurw 



m— 1 



2 y /N 1 N 2 J (e 2iXcos ^ r+ ^ - l) cos (ifo - ifo + Qr)«i(r, t)u 2 (r, t)dr 

m—oo , . \ \2m+l r 



V + (e^-^ + e-W^-v)) = 2z VA^V2 Ji(2A) cos - ^ - p). (60) 

z — '„ (2m + ll! mlim + 1)! V / 



(2m + 1)! m!(m+ 1) 



711 = 



Substituting ([55]) and ([6H)l into ([ST)) , and doing the integral over (/?, we finally obtain the central result of this 
section: 



Z(iA;%/iVie #1 ,%/iV 2 e lV;2 ) = / * ^ e (Jo(2A)-l)(M+JV2)+2iV77IW?Ji(2A) cos _ (g-n 

e (Jo(2A)-i)(Jv 1+ Jv 2 ) Jq ^2 A /iViiV2J 1 (2A)) . (62) 



B. Interference of independent clouds in number states 

Let us now explain how to calculate Z' (iX, Nx, N 2 ) for the Fock states with the number of particles equal to Nx 
and N 2 . First, we need to e xpand the Fock states \Nx, N 2 ) using the coherent states basis. Since the basis of coherent 
states is overcomplete [l08l |. there are many ways to do a decomposition. For our purposes it is convenient to use 

\Nx,N 2 ) = (a fl!l 2 10) = ^/NxANj.a-^-^e* 2 l** d ^ Rd ^ 2R e ~^ iyiR -iiv 2y2fl | Qe » yiR) ae W (63) 

vW v%! Jo Jo ( 27r ) 

where a is an arbitrary real positive number. Coherent states are given by 

\ae itpiR ,ae ltp2R ) = e -<* > +<*e ivl *4+<*e iv ™<4\o) i (64) 
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and the overlap between them equals 

(ae tlpiL ,ae ltp2L \ae lipiR ,ae llfi2R ) = ^{^R-^+^r-^l-^ _ ^ 

One can also expand the bra state {Ni , N 2 | similarly to (|63)l by introducing integration variables p.\L , P2L ■ For any 
given values of tp, <Pil, <^2L, Pir, and tp 2 R, matrix elements between coherent states can be evaluated as in previous 
section, and by simple modification of Eq. (|61|) we obtain 



zf(ix,N U N 2 ) = J J J J j*e*£ii*p^ 

(ae lVlL ,ae lip2L \ae lipiR ,ae tip2R ) x 

One can now substitute Eq. (|65p into the equation above, introduce variables 

<pi = Vil - Pir and ip 2 = ifi2L - P2R, (67) 
and integrate over p, tpm, P2R- Multiple cancellations occur, and eventually we obtain 

Zf(iX,N u N 2 ) = J J ^^JV 1 !JV 2 !a- a ^- 2 ^e Wl ^^ a ^e Jo ( 2A )° a ( e " i * ,1 + e " iw )Jo (2a 2 J 1 (2A) e - < ^ 1 +^)/ 2 ) .(68) 

Both integrations in the equation above can be done in a closed form for arbitrary positive integer N\ and N2 using 
hypergeometric functions. Here we will present the results only for N\ — N 2 = N. One needs to expand the last 
exponent and Jo (2a 2 Ji(2X)e~ 1 ( ipi+V2 ^ 2 ) using Taylor series. After integration over dpi and d<p 2 dependence on a 
disappears, and we obtain the final result for Fock states: 

* *) = tm i Jo % )n (iJl ^^~ k) = 2 * ~ n > 1; ) Mm™ > ^ 

where 2F1 (a, 6; c; x) in a hypergeometric function defined by 

7-1/ 7 s 1 ■ Q&!C , o(o + l)6(6+l) a; 2 
2Fl (<2 ' 6; 05 = 1 + T II + c(c+l) 2! + - (?0) 



C. Clouds with a well denned relative phase 

Let us now consider the case of clouds with a fixed relative phase, when the initial state \(fio,N) at t — is given 
by Eq. ([15]). This state can be expanded using the coherent states basis as 

l^o,A0 = lnM * „ (ole^/ 2 + 4e- <yo/2 ) Jv 10) - VW.(V2a)~ N ^- e ~ lNtfiR \ae^ R+ ^^ 2 ,ae lVR - lv °/ 2 ). 
(2" NY) 1 " Jo 2n 

A similar expansion can be written for bra- vector (po,N\ using the phase variable tpL- Coherent states and their 
overlaps are given by Eqs. (|64|) -(|65 p . and one obtains an expression for the generating function Z(iX, N) similar to 
Eq. ((Ml): 

Z(iX,N) = J J J ^^ A^ 



Xf 



2{J {2\)-l)a 2 e iv R- iv L io? J x (2 A) (e i(v " ™ +f R ~v L ) +e -i(<f-Vo+VL -<pr) ) 



The integrand depends only on the difference Ap — pl — (pR, and the integral over <p can be done analytically. 
Dependence of Z(iX,N) on po drops out, as expected: 

Z(iX, N) = J J ^Ma- 2N 2- N e mA *e 2J °^ a2 ^J (2a 2 J 1 (2X) e -^) . (71) 

Expanding the last exponent and Jo (2a 2 Ji(2X)e~' tA ' p ) in the expression above and integrating over Ap, we obtain 
the final expression for even N : 

z<a, W) . W!2 -» |««^ " 2 - {M) • < 72 > 

where C"(x) is a Gegenbauer polynomial [lO 1 
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FIG. 8: Experimental setup for interference of 2D gases. Note that the interference patterns are straight at low temperatures 
indicating suppressed phase fluctuations. Meandering patterns at high temperatures come from strong phase fluctuations. 
Taken from Ref. 1321. 



IV. INTERFERENCE OF FLUCTUATING LOW-DIMENSIONAL GASES 



As we discussed in previous sections, for macroscopic three dimensional Bose-Einstein Condensates the long-range 
phase coherence manifests itself in the nearly perfect interference fringes between two independent condensates [27| . 
For low dimensional Bose gases, the situation is different, since phase fluctuations are very effective in destroying the 
long-range order. In one dimension, long-range coherence is prohibited even at zero temperature 92], while in two 
dimensions any nonzero temperature destroys long-range order [9lj ]. In addition, the Berezinskii-Kosterlitz-Thouless 
(BKT) phase transition occurs [85l |86|. which separates the low temperature phase with power-law correlations from 
the high temperature phase with short-range correlations. Phase fluctuations reduce the average visibility of the 
interference fringes, and result in the shot to shot fluctuations of the visibility. 

In this section, we discuss how measurements of interference fringes can reveal information about spatial correlations 
within individual condensates. The typical experimental setups are shown in Figs. 2] and [5] They correspond to the 
so-called open boundary conditions (OBC). Essentially the OBC mean that the imaged area is cut out of a larger 
system. As a theoretical model one can also consider a one-dimensional condensate with periodic boundary conditions 
(PBC), which corresponds to interference experiments with two coaxial rings lying in two p arallel xy planes. While 
this model is somewhat artificial from the experimental point of view (see however Ref. [l09j |). it allows a very elegant 
theoretical analysis, hence we will discuss it in these lecture notes as well. 

The confining potential is highly anisotropic, and after it is switched off, the clouds predominantly expand in 
the transverse direction, while no significant expansion occurs in the axial (for ID gases) or in-plane (for 2D gases) 
directions. For low-dimensional gases the phase of the condensate doesn't have a long-range order due to quantum 
or thermal fluctuations. Locally the phase determines the positions of the minima of the absorption intensity, and 
fluctuations of the phase lead to fluctuations of the interference fringe positions along the condensates, as shown in 
Fig. [5] Fluctuations of the fringe positions contain information about the original phase fluctuations present in the 
system, which are preserved during expansion. 

To extract information about fringe position fluctuations for the ID case, we will integrate the intensity along the 
axes of the clouds. Fluctuations of the relative phase result in fluctuations of the minima positions for different x. 
For each y, the image can be integrated along the x direction to give the integrated fringe amplitude Pq{L) (see Fig. 
[4|. Note that the integrated fringe amplitude depends on the integration length L. One experimental image can be 
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used to extract information for different values of L. Many images are still required to obtain distribution functions 
for each L. For 2D gases, the setup is analogous and is shown in Fig[S] Here part of the integration is performed by 
the imaging beam itself. The size of the integration area along the direction of the imaging beam can be controlled 
by applying magnetic field gradients, so that only a specified section of the cloud is resonant with the probe light. 

The operator which corresponds to the fringe amplitude pq(L), illustrated in Fig. |4l is the same as pq defined by 
Eq. |27|) . where the integration along x— dimension is limited to the section of length L. Let us first consider the 
expectation value {\pq{L)\ 2 ) (expectation values of operators which depend on the phase of Pq(L) vanish, similar 
to 3D case, since two clouds are assumed to be independent). One has to use modified formulas (|3T)|) - ([3"Tj) , where 
operators a\, a,i are now allowed to have x— dependence. In the limit when the number of particles in the section of 
size L is large, the average value of (\pq(L)\ 2 ) is given by [5f| 

(\PQ(L)\ 2 )= / dxxdx'Mixx^Ux'Jaxix'Jaiixt)). (73) 
Jo Jo 

Note that in Eq. (f73|) we used the normal ordered form of the operators which means that we neglect the shot 
noise considered in section IIIII This is justified for long condensates as we discuss below. ^From now on, we will 
concentrate on the case when independent clouds are identical and have the same density of particles with equal 
interaction strengths. Then 

(\PQ(L)\ 2 ) = [ L [ L dx 1 dx' 1 (a^x 1 )a(x' 1 )) 2 . (74) 
Jo Jo 

To gain intuition into the physical meaning of the average amplitude of interference fringes, we address two limiting 
cases. First, consider the situation when (aJ (x)a(0)) decays exponentially with distance and the correlation length 
is given by £ << L. Then Eq. (|74|) implies that |/Oq(L)| oc y/L^, which has a simple physical interpretation. Since 
the phase is only coherent over a length £, the system is effectively equivalent to a series of L/£ pairs of independent 
condensates. Each pair contributes interference fringes with a constant amplitude proportional to £ and a random 
phase. The total amplitude pq (L) is the result of adding L /£ independent vectors of constant length £ and random 
direction. Adding random uncorrelated vectors gives a zero average except for a typical square root fluctuation. Thus 
scaling of the absolute value of the net interference amplitude is y/L£. This observation is similar in spirit to that 
made in Ref. [65j of interference between 30 independent condensates in a chain. Fringes can be seen, though their 
average amplitude is suppressed by a factor of v30 compared to the interference between two condensates. Now 
consider the opposite limit of perfect condensates, for which (a> (x)a(0)) is constant. In this case Eq. (|74f implies 
that |pq(L)| oc L. Pictorially this is the result of adding vectors which are all aligned, resulting in a fringe amplitude 
absolute value of which scales as the total size of the system. 

Methods developed in this section for analyzing \pq(L)\ 2 can be applied to condensates with either uniform and 
non uniform densities. For simplicity, we concentrate on the case when L is much smaller than the size of the clouds, 
so the change in the atomic density along t he clouds can be ignored. In this case correlation functions for ID gases 
are described by the Luttinger liquid theory jllOUlllI ]. For OBC at zero temperature two-point correlation functions 
are given by 

(al(x)a(y))~p(Z h /\x-y\f 2K . (75) 

Here p is the particle density, £/j is the healing length, which also serves as the short range cutoff, and K is the so- 
called Luttinger parameter, which characterizes the strength of interactions. For bosons with a repulsive short-range 
potential, K ranges between 1 and oc, with K — 1 corresponding to strong interactions, or "impenetrable" bosons, 
while K_—*oo for weakly interacting bosons. Substituting Eq. (f75|) into Eq. (|74|) and assuming that L > ^, we 
obtain [55[ 

c \ 1 / K 



(\p Q (L)\ 2 )=Cp 2 L 2 ^f-) , (76) 

where C is a constant of order unity. We see that the amplitude of the interference fringes (\pq(L)\) scales as a non 
trivial power of the imaging length. In the non interacting limit {K — ► oo) the scaling is linear (\pq(L)\) ~ L, as 
expected for a fully coherent system. Interestingly, (\pq(L)\) ~ \J~L appears in the hard core limit (K — 1), as in 
systems with short range correlations which were discussed above. 

One may be concerned that Eq. (|76p gives only the long distance asymptotic behavior of the correlation functions, 
and does not describe the short distance behavior. From Eq. ([74]) one finds that the contribution of the short 
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distance part of the correlation functions to (\pq(L)\ 2 ) scales as L. In the physically relevant case of K > 1 and in 
the limit of large L this contribution is smaller than Eq. ([75]) . We note that in principle one can use t he ex act Bethe 
ansatz solution of the Lieb-Liniger model |ll2| to obtain correlation functions valid at all distances [l!3| . Another 
contribution which has been neglected is the shot noise. The shot noise contribution to (\pq{L)\ 2 ) comes from the 
normal ordering of operators, and equals 

dxidx' 1 {a\{xi)a2{xi)a\{x' 1 )ai(x' 1 )) / dxidx' 1 (a\(xi)a\(x' 1 )ai(x' 1 )a2(xi)) = 



o Jo Jo 

L r L 



o Jo 



dxidx' 1 S(xi — x' 1 )(a\(xi)ai(x / 1 )) = ni£>L. (77) 



In the limit of large L and K > 1 this is again a subleading contribution and can be neglected. 

For 2D, one can use similar approach to describe the contrast distribution at finite temperature below the BKT 
transiti on. We n ote that we assume that the temperature is small enough such that 2D gas is in a quasiconde nsate 
regime |ll4Ull5j . when only phase fluctuations are present. In this case, correlation functions are given by (85l. [86l . fll6j 

(at(r)a(0)>~p(^J , (78) 

where rj(T) = mT / '(2irh 2 p s (T)) depends on the temperature and the superfluid density p s (T). The BKT transition 
takes place at the universal value rj c (T c ) — 1/4. To keep connection to the ID case, we introduce 

K = 1/(277(T)), (79) 



and restrict our attention to K > K c = 2. For temperatures above the BKT transition Eq. (|78|) doesn't hold, and 
correlations decay exponentially. This means that the integrated interference amplitude will only increase as the 



square root of the integration area 55 1 



Fig. [8] illustrates experiments performed with 2D gases to identify the BKT transition [32|. Two independent 
2D condensates are confined in transverse directions using an optical lattice potential. After the optical potential is 
switched off and clouds expand, the density is imaged on a CCD camera. When temperature is small, interference 
fringes are straight lines. As the temperature is increased, the fringes start to meander due to spatial fluctuations of 
the phase. Integrating the image along the section of variable length L in x— direction gives L— dependent fringe 
amplitude {\pq{L)\). Scaling of this amplitude with L contains information about rj{T), which is expected to have a 
universal value i] c (T c ) — 1/4 at the BKT transition. 

Fig. [9] illustrates the procedure used to extract scaling exponents in experiments. (\po(L)\ 2 ) (denoted as (C 2 ) in 
Fig. [9K) is plotted as a function of L, and its scaling with L in a certain range (see Ref. [32| for more details) is used 
to extract the exponent r/(T) (denoted as a in Fig. 03). The average central contrast Co serves as a "thermometer", 
such that smaller values of Co correspond to higher temperatures. Above the BKT temperature the value of r](T) 
extracted from interference experiments is expected [55j to be equal to 0.5, while at the transition point it is equal to 
Vc (T c ) = 1 /4. Fig. [9)3 shows a sudden change of the exponent in a relatively narrow rang e of temperatures. This change 
is reminiscent of the universal jump in the superfluid density for 2D helium films [117j |. Of course, since experiments 
are done with finite systems and imaging done along the y— direction performs averaging over the inhomogencous 
density profile, one shouldn't expe ct the uni versal jump, but rather a crossover. The presence of the trap also affects 
parameters of the BKT transition |l!8l . Ill9| . 

The BKT transition is an example of a topological transition which is driven by the unbinding of vortices, and the 
remarkable feature of experiments in Ref. [321 ] is the ability to independently resolve the vortices. When only one 
of the condensates has a vortex, the interference pattern will have a disclination-like structure [31( . It was shown 
experimentally in Ref. [32j that proliferation of vortices occurs at the same point in the parameter space as the jump 
in the scaling exponent. In following sections we will show that in uniform systems not only the scaling exponent, 
but also the full distribution function of the fringe contrast has a universal form at the BKT transition. 



A. Interference amplitudes: from high moments to full distribution functions 

Measuring atom density to obtain an interference pattern is a classical measurement on a quantum mechanical wave 
function. The process of the measurement itself introduces an intrinsic quantum mechanical noise. Said differently, 
from shot to shot we will not get precisely the same value of the amplitude of interference fringes. Expressions for 
(\pq(L)\ 2 ) which we derived in Eqs. (|73p ~ (l76p correspond to averaging over many shots. For example, data points 
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FIG. 9: Emergence of quasi-long-range order in a 2D gas. a, Examples of average integrated interference contrasts (C 2 ) are 
shown for low (blue circles) and high (red squares) temperatures; L x is the integration length along x— direction. The lines 
are fits to the data using the power-law functions l/(L x ) 2a . b, Exponent a as a function of the central contrast Co. Central 
contrast Co serves as a "thermometer", such that smaller values of Co correspond to higher temperatures. Dashed lines indicate 
theoretically expected values of a above and below the BKT transition in a uniform system. Taken from Ref. [32] ]. 



from Fig. [9j^ correspond to averaging over approximately one hundred measurements (33j . However each individual 
shot will give the value of |pq(L)| 2 which may be different from its average value. 

An interesting question to consider is how this amplitude fluctuates from one experimental run to another. To 
address this question we need to consider higher moments of the operator |pq(i)| 2 . Generalizing the argument which 
lead to Eq. ([74)) , we obtain 

(\PQ( L )\ 2n ) = ••■/ dx\ . . . dx n dx[ . . . dx' n \(a) (xi) . . . a t (a;„)a(x' 1 ) . . . a(x' n ))\ 2 . (80) 
Jo Jo 

In Eq. ([BH]) we used a normal ordered correlation function similar to Eq. ([74"]) . One can calculate [99| corrections 
due to normal ordering for higher moments of |pq(L)| 2 , and show that in the limit of large L and K > 1 they can be 
neglected. 

^From Eq. (|74j) we observe that (|/?q(£)| 2 ) contains information about two-point correlation functions of individual 
clouds. Eq. ([50)) shows that higher moments of |/9q(£)| 2 contain information about higher order correlation functions. 
The full distribution of the fluctuating variable |pq (L) | 2 contains information about all high order correlation functions. 

In the Luttinger liquid theory fluctuations of the phase are described by the Gaussian action. For G auss ian actions 
higher order correlation functions are simply related to two-point correlation functions (see e.g. Ref. [l!6j |V. 

UiaHxMx'j)) 

(atfo) . . . at(^)aK) . . . a(z'J) = n ^ -gL_ )} ■ (81) 

i<j i<j 

Using this formula together with Eq. (|75[) . the higher moments of fringe amplitudes can be written as 



(\pQ(L)\ 2n ) = A 2 n Z 2n , where A = Jc P 2 ^ h /K L^/ K , (82) 
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C is a constant of the order of unity and for OBC in ID 



Z, 



dui...dv n 



U.\ u i ~ Uj\H\vi - Vj 

i<j i<j 



U\ n i- v j\ 



Here for OBC G(x, y) is given by 



o 



G(x,y) = log\x - y\. 



du!...dv n e \ t<3 l<3 ,J /.(83) 



(84) 



Integrals similar to Eq. (|8"3"|) appeared in the literature before tl20L 121 1 , b ut they are not easy to compute. 

The Gaussian model possesses a powerful conformal symmetry llll . 116j | , which dictates the form of the correlation 
functions for periodic boundary conditions or nonzero temperatures. For PB C wit h circumference of the condensates 
equal to the imaging length, the change in the correlation functions leads to 111] 



G per (x, y) — log — sin7r|x — y\. 



(85) 



For nonzero temperature, Z 2n depends on K and the thermal length £t = Tw s / {ksT), where v s is the sound velocity: 

G(x,,,|).log(gsinh^-^). (86) 



The analysis in this section (except for section llV C| doesn't depend on the particular form of G(x,y). The only 
general restriction is 



G(x,y) = G(y,x). 



(87) 



Eq. (|83j) is also valid for 2D case below the BKT transition temperature, with m and Vi being 2D variables on a 
rectangle with G(x,y) = log \x — y\, and with a properly redefined A a . 

In what follows, we will be interested in distribution functions PF(a) of the variable a = \pq(L)\ 2 /Aq and of its 
normalized version a — \pq(L)\ 2 / (pq(L) 2 } . By the definition of higher moments 



W(a)a n da, and — ^ 



W{a)a n da. 



The general problem we consider now is how to construct distribution functions W(a) given by Eq. (|88p with 
Z 2n given by Eq. (|83[) . G(x,y) can be an arbitrary symmetric function of x and y, and not necessarily the function 
of their difference. This allows us to stud y systems with non-uniform density in external traps, where one needs 
to use modified correlation functions [l!5j . In section HV Bl we will show that W(a) is connected to the partition 
functions of various Sine-Gordon models and Coulomb gases, and methods developed in this section provide a new 
non perturbative tool for calculating partition functions of such models. 

One can think of Z2 n in Eq. ([55)) as partition functions of a classical two-component gas of fictitious charged 
"particles" in a microcanonical ensemble, with K being the "temperature". We first briefly comment on the limiting 
cases K 3> 1, and K — > 1. For K ^> 1, expansion of Z 2n gives Z 2n ~ 1 in the zeroth order. This means that 
W(a) w S(a — 1), i.e. a very narrow function peaked near its average value. Higher order terms in the expansion give 
W(a) a small width of the order of 1/K, and are studied in detail in Appendix [SJ For ID gases at zero temperature 
with K — > 1 or nonzero temperatures with £tK/L <C 1 the distribution function W{a) is Poissonian irrespective of 
OBC or PBC: W^(S) = e~ a . To demonstrate this we need to verify that Zi n fZ 2 = a n e~ a da = nl. This can be 
shown using the classical gas analogy: as K — > 1 , Z 2 — J l/|rc — y\^' K dxdy starts to diverge for x — > y, and the main 
contribution to Z2 n comes from "molecular" states of the two-component gas, i.e. from the parts of the configuration 
space in which each "particle" has a "particle" of the opposite charge in its neighborhood; n\ is just the number of 
ways to form such pairs. In this language K 3> 1 corresponds to the "plasma" phase of the classical charged gas, 
and the evolution of the distribution function can be understood as a formation of "molecules" out of the "plasma" 
phase as the "temperature" (i.e. the Luttinger parameter K) is lowered. For finite £t, the main contribution to Zi n 
comes from distances ~ ^tK/L, and if this parameter is much smaller than 1, then again "molecular" contributions 
dominate. 

There is a simple physical interpretation to W(a) for both K 3> 1 and K — > 1. When bosons do not interact and 
K = oo, there should be no phase fluctuations within individual condensates. Hence in each experiment we should 
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find a perfect interference pattern although positions of the density minima are unpredictable (see discussion in section 
HI)) . Alternatively in the regime of strong interactions when K — > 1 we can think of the net interference as a result of 
adding many random uncorrelated two dimensional vectors (see discussion after Eq. ([74]) ). Earlier we used the fact 
that for a random walk the net displacement is proportional to the square root of the number of steps. But we also 
know that for 2D random walks the distribution function of the square of the net displacement is Poissonian, which 
is what we find for W(a). 

As K varies from 1 to oo, the distribution function W(a) should evolve from being a very broad Poissonian function 
to a narrow delta function. The evolution of the distribution in the intermediate regimes will be studied in detail 
below. In Appendix [A] we develop a systematic expansion of Zi n in powers of 1/K for large K which works for different 
dimensions and boundary conditions, and investigate the distribution functions in this limit. The rest of this section 
is organized as follows. In section [IV Bl we discuss the connection of our probl em t o the Sine-Gordo n mo dels, which 
describe various physical problems, ranging from flux line in superconductors jl22| to string theory [l23{ . In section 
IIV CI we obtain full distributions for the ID case with PBC using the exact solution of the boundary Sine-Gordon 
model on a circle [56[ . In section IIVDI we present a novel non perturbative solution which is applicable for any value 
of K and for various dimensions and boundary conditions. We discuss the connection of the distribution functions 
of fringe visibilities to the statistics of random surfaces 58], and prove that for ID case with periodic boundary 
conditions in the limit of large K the distri bution of fringe visibilities is given by the Gumbel distribution, one of the 
extreme value statistical distributions Il24ll . 



B. Connection of the fringe visibility distribution functions to the partition functions of Sine-Gordon models 

To illustrate the connection of the fringe visibility distribution functions to the partition functions of Sine-Gordon 
models we start from a ID system with periodic boundary conditions. Shortly we will demonstrate that this example 
provides a unique possibility to make an excursion into very different subjects of physics and mathematics. As shown 
earlier, for a one-dimensional ring condensate with PBC and circumference equal to the imaging length, the correlation 
functions lead to 

G per (x, y) — log | — sin(7r(a; — y))\- (89) 

This function describes a two-component gas of particles interacting via a ID Coulomb potential, which is a logarithmic 
in the interparticle separation. For PBC the distance is a chord function sin(7r|a; — y\)/ir. It is common to define the 
grand partition function of such Coulomb gas as 

00 2n 

z °(^) = E£n2^r> (9°) 

where is given by Eq. with G(x,y) substituted by G per (x,y). 

Several different physical problems can be related to t he p artition function given by Eq. (|90p . Here we mention 
only a few examples: a) the anisotropic Kon do m odel [l25j : b) the quantum impurity model in one dimensional 



interacting electron system s int roduced in Ref. [lOOf , which has been extensively studied in the context of edge states 
in Quantum Hall systems [l26j ] ; c) the b ackgr ound-indep ende nt string theory and the model of strings attached to 
a D- brane originally introduced in Ref. jl23l | (see Ref. jl27| for a recent review); d) Calogero- Suthe rland model 
128], which has numerous application as an effective model; e) flux line pinni ng in superconductors [l22j |: f) quantum 
tunneling in the presence of dissipation within Caldeira-Leggett model [129| ; g) interference of two one-dimensional 
condensates [55l l56l. [58j. 

We now briefly describe how the expression given by Eq. (|90[) appears as a partition function in physical systems. 
We consider the imaginary time action 

S P er[g] = \ I dx f dT[{d T <f) 2 + {d x <j)f] + 2g f dTCOs[l3(j>{x = 0, t) + 2tt P t] (91) 

^ J-oc JO JO 

known as the boundary Sine-Gordon model [l30j |. The quantum field 4>{x,t) in Eq. (|9ip is defined on an infinite 
line in x direction and is assumed to be periodic along r direction: <f>(x, r) = 4>(x,t + 1). The interaction term is 
present only at x = 0. A typical physical system which is described by action (f9Tj) is an interacting ID electron liquid 
scattered by an impurity [lOfJ. In this case the cos-term describes backscattering of electrons within the Luttinger 
liquid formalism. In Eq. (|91[) we also added a p— depe nden t phase winding term. For a quantum impurity problem 
this is somewhat reminiscent of having a finite voltage jl3l| , while in interference experiments this term corresponds 
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to a relative momentum along the condens ates (56j (except for section IIV CI we only consider the case with p = 0). 
Partition function of Eq. (f9"Tj) is defined as |l20l ] 

One can expand Z p (K,g) in Taylor series of the coupling g. Nonvanishing contributions come from combinations 
which have equal number of exp(+i/3<^)) and exp(— if3(f>) terms. This is essentially a charge neutrality condition for the 
Coulomb gas. If we identify 

3 2 1 

^=K> x = g(2n)^, (93) 

we obtain 

oo 2n 

W5) = EW^ ) - (94) 
Here the microcanonical partition functions z!$ are represented by the following integrals 



2n ~Jo "Jo l\2n2n e 
Zifo for p = is related to Z^ defined above by 

= Z^(2n)~^ K . (96) 

Then 

°° 2n 00 In °° 2n r<x> poo 

^^Et^^Et^KHEt^ / W(a)a n da = W{a) I {2g^) da. (97) 

n=0 ^ '' n=0 * '' n=0 ^ '' J J ° 

Zo(K,g) is essentially a Hankel transformation of W{a), and inverting Eq. (|97p we can express the probability 
W(a) through the partition function Z (K,g). Noting that Io(ix) = Jq{x) and using the completeness relation for 
Bessel functions, J Q °° Jq(Xx) Jo(Xy)\x\ XdX — S(\x\ — \y\), we obtain 



n . < . 2s in(^)n fe<i 2sh 



a, fe 2sin (^) 



(95) 



W(a) = 2 / Z Q (K,tg)J (2g^)gdg. 
Jo 



It is important that the last equation has the partition function at imaginary value of the coupling constant. This 
should be understood as an analytic continuation of Zo(K,g). For PBC in ID, partition functions Z p (K,g) can be 
evaluated using the exact solution of the boundary Sine-Gordon model with periodic boundary conditions, and this 
approach will be presented in section IIV CI 

For open boundary conditions one can also write the grand canonical partition function 

00 2n 

^) = E^». (99) 

n=0 * 

where Z 2n are given by Eqs. ([53")) , We can express Eq. as a partition function of a certain Sine-Gordon 

model: 

where 

/OO p OO p 1 

dx j dr[{d T (p) 2 + {d x <j)f] + 2g / dr cos[27T0(ie = 0,r)]. (101) 
-oo J — oo J 



24 



Note different limits of the r integration in the first and the second terms. One can see that Eqs. (|99| and (|100|) 
define the same Z(g) by observing that — log (\x — y\) / (2tt) is a free propagator of the Gaussian action on a plane. 

So far we established the relation between the distribution functions of the amplitude of interference fringes of ID 
condensates and the boundary Sine-Gordon models (|9"Tj) . (|10ip . Generalization of this argument to the case of 2D 
condensates is straightforward. Higher moments of the interference amplitude are given by integrals of type (|83[) . but 
each u and v is now a two-dimensional coordinate. Distribution function of fringe amplitudes is then related to the 
partition function of the bulk Sine-Gordon model 

/OO />OG />1 />1 

dx dT[{d T <t>) 2 + {d x <j>) 2 } + 2g / / dTdxco8[2ir^(x,T)]. (102) 
-oo J-oo JO JO 

Expanding the partition function corresponding to action (1102)) in powers of g, one finds the expression identical to 
Eq. (|97|) . but with ty(a) corresponding to interference of 2D condensates. 

When we describe interference of systems with OBC, we use Eqs. (llOip and (|102p . in which the field 4>(x,t) is 
defined on a whole plane, and is not periodic in r. In both these cases the interaction is present only in some part 
of the system, so translational invariance is lost in both x and r dimensions. Thus to calculate using Eq. 

(|98|) one needs to calculate partition functions of inhomogeneous Sine-Gordon models. Exact solution of Sine-Gordon 
model is available only for periodic boundary conditions, so to treat open boundary conditions one needs to develop 
alternative methods. Conversely, if one has a solution for W(a), this provides a tool for calculating partition functions 
of inhomogeneous Sine-Gordon models and Coulomb gases using Eq. I|97p. In scction UVDI we present a novel mapping 
of W(a) for arbitrary G(x, y) to the statistical properties of random surfaces, which provides a new tool for calculating 
partition functions of a wide class of Sine-Gordon models and Coulomb gases. In general, partition functions with 
inhomogeneous g[x, r) can be evaluated using this mapping as well. In this case the function g(x,r) will appear 
in the integrand of Eq. (|102p . We point out that the mapping which we use is not related to the existence of the 
exact solution of Sine-Gordon models, but relies only on the structure of the correlation functions in the absence of 
interactions. We also note that a suitable extension of our method can be used to compute correlation functions of 
Coulomb-gas models in equilibrium and non-equilibrium situations. 

C. Distribution functions for ID gas with periodic boundary conditions 

There are three natural ways to compute partition function Z$(K, g) for complex g, which can be used to construct 
distribution functions for ID gas with PBC using Eq. ((98]) . The first one is related to the theory of symmetric 
polynomials and is described in Appendix [Bj The second one relies on the integrability of the quantum impurity 
model defined by the action (|9ip . and makes use of the thermodynamic Bethe ans at z (T BA) . This approach is 
discussed in Appendix [Cl These two approaches have been presented in Refs. [l20j |. [l32]. They are difficult to 
implement directly and are included for completeness. The third approach described below in section II V C II leads to 
the most transparent answer. Although it is intrinsically related to the first two, it has a broader applicability and 
will be studied in detail. 

At imaginary g the theory (|9ip is apparently non-Hermitian, but as we will see later the model belongs to a special 
class of non-Hermitian field theories, which have real spectrum. We describe an interesting connection with VT- 
symmetric quantum mechanics and conformal field theo ries with negative central charges. Here we also make close 
contact with the impurity problem in a Luttinger liquid [lOOj . This connection is based on the following observation: 
when we consider Taylor expansion of the quantum impurity partition function in powers of <?, we need to take multi- 
point correlation functions at the same x but at different times. To calculate higher moments of the interference 
amplitude signal we need to use correlation functions for equal tim e, but at different points in space. The two are the 
same because of the relativistic invariance of the Luttinger liquid jllll . Illd ] . 

1. Mapping to integrable structure of CFT and singular anharmonic oscillator 

To introduce our method we need to make several formal remarks regarding one dimensional conformal field theories. 
In a recent series of papers Bazhanov, Lukyanov, and Zamolodchikov expl ored an integrable structure of conformal 
field theories focusing on connections to solvable problems on lattices [l33l |. The key ingredients of the solvability of 
lattice models are the so-called transfer matrix operators T(A). These operators contain information about all integrals 
of motion as well as excitation spectra of the system. Transfer matrices are defined as a function of the so-called 
spectral parameter A (in the continuum limit A corresponds to rapidity) and commute for different values of A. The 
latter property is a direct manifestation of the existence of an infinite number of commuting integrals of motion. In 
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his studies of the 8- vertex model, Baxter jl34j introduced the operator Q(A) which helps to find eigenvalues of T. 
Operators T and Q± satisfy a set of commutation relations, in particular [1331 ] 

T(A)Q±(A) = Q±(qX) + Q ± (q- 1 X), (103) 

where q — exp(iir/2K). So T matrices can be obtained explicitly when one knows the Q operators. 

By construction, these operators act in the representation space of Virasoro algebra which can be constructed from 
the Fock space of free bosonic operators a± n defined such that 

ft 27T 

a n \p) = 0, for n > 0, P\p)=p\p), [a n ,a m ] = -—5 n + m ,o- (104) 

z A 

Here \p) denotes the vacuum vector and p is the zero mode of the corresponding bosonic field 4>{u): 

0(«) =iQ + iP u + J2—e mu . (105) 

Operators A ± (A) = Q ± (A)A=f 4 p/^ act in the representation space of Virasoro algebra, which can be constructed from 
the Fock space of bosonic operators a± n satisfying a n \p) =0, (n > 0). The Fock vacuum state \p) is an eigen state of 
the momentum operator, P\p) = p\p). For p = N/2 (N = 0, 1, 2...) the vacuum eigenvalues of the operator A±(A), 
A± (A) \p) = A± ac ^ \p) , are given by (below we consider only the quantities with the + subscript which correspond to 
the positive p) 

A( vac 1(\) = Z p (K,-ig), (106) 

where 

*-jf-(£). (107) 

At this point it is not clear what we gained by connecting the analytically continued partition function of the impurity 
problem Z p (K,ig) to the expectation value of the operator A^ vac \X). As we discuss below, a considerable number 
of important results have been derived for A^ vac \X). We will be able to make use of these resul ts to obtain the 
distribution functions of the fringe amplitudes. The function A^ vacS) has known large-A asymptotics |l33| 

log^^XA) ~ M{K) (-A 2 ) , (108) 

where the constant M(K) is given by 

M(K) = ^ { ^-i m r( 2K K \ ■ (109) 

The function J 4^ oc - ) (A) is entire function for K > 1 and is completely determined by its zeros A&, k = 0, 1, .... Therefore 
A<- vac ^{\) can be represented by the convergent product 



00 / \2\ 

A {vac) (X) = T] 1-T2 ' ^ ( " ac) (0) = 1. 



(110) 



On the basis of analysis of a certain class of e xactl y solvable model, corresponding to the integrable perturbation of 
the conformal field theory, it was conjectured in [l35j | that the so-called F-system and related T system (where Y — e £r 
and e r are the Bethe-ansatz energies parametrized by r, the nodes of the Dynkin diagrams (see TBA-section) satisfy 
the same functional equations and possess the same analytical structure and asymptotics as the spectral determinant 
of the one-dimensional anharmonic oscillator. Further, the same functional equations, analytical properties (|110p and 
asymptotics (I108[) are satisfied for the vacuum eigenvalues of Q-operator for special values of p and the latter are 
given by the spectral determinant of the following Schrodinger equation 

(-dl + x 2a )t>(x) =EV(x). (Ill) 

The spectral determinant is defined as 
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Soon after, in Ref. 136], this conjecture has been extended to all values of p : 

A^(\,p)=D{p\ 2 ), (113) 
where now D(E) is the spectral determinant of the Schrodinger equation 



- %9(x) + [x iK - 2 + ^r^J *0c) = E*{x), (114) 

with I = 4pK—l/2, Here p = (AK) 2 ~ 1 / K [T(1-1/(2K))] 2 . In relation to the interference problem, p = (md/ht) tan 9, 
where m is the atom's mass, d is the separation between the two condensates and t is the time when the measurement 
was done after the free expansion started. 

For some values of parameters the Eq. (|114[) can be solved exactly: 

a) K = 1. In this case, corresponding to a singular harmonic oscillator E n = 4n + 21 — 1, n = 1, 2, ... and 



r(| + i)e- GE 



D(E,l)= Z\ 2 ' E . , (115) 



. 4 1 2 4 > 

where C is nonuniversal renormalization constant. 

b) The case if — > oo is recovered by the rigid well potential for which the eigen energies are given by the zeroes of 
the Bessel function, and therefore 

D(E, I) = T(l + 3/2)(Ve/2)- 1 -^ 2 J 1+1/2 (VE). (116) 



c) For K — 3/4 and I — the result is expressed in terms of Airy function. 

For generic values of if and I the Schrodinger equation can be solved numerically with a very good precision with 
subsequent computation of the spectral determinant. Alternatively, fo r n > 5 — 10 the spectrum of the equation (|1 14[) 
is very well approximated by the standard WKB expression (see e.g. [l37l |) 

E n = e(K){n- 1 i(K))^, (117) 

where n = 1,2, ... Here 71 (if) is the Maslov index. For 1/2 < K < 00, ji(K) = | - |, for K = 00, 7* (If) = -1/2 and 
for < K < i, 7/ (if) = 4K gg~ l . The function e(if) in Eq. (117]) reads 

< K ^= r(i + ^) 

In principle, the function 7; (if) can be a smooth function interpolating between limiting values given above and 
can be considered as a noninteger Maslov index jl38l |. This interpolation allows to (approximately) evaluate the 
partition function and the distribution function in many cases. In the limiting cases K — > 1 and K — > 00 the WKB 
approximation gives the exact spectrum. We point out however, that using the approximate WKB spectrum carelessly 
can result in non-physical results, e.g. negative values of the distribution function W(a), and thus have to be used 
with care. WE emphasize that the solution of the ODE (I114j) for several potentials gives us already analytically 
continued function Z p (K, ig). It seems that this approach to finding the solution is easier and more elegant than the 
solution of TBA equations with subsequent analytical continuation. Moreover a variety of approximate methods are 
available for solving Schrodinger-like ID equations. 

Note that Eq. pi4[) has an interesting duality symmetry which generalizes the Coulomb-harmonic oscillator duality 
and which allows to relate the K > 1/2 and < K < 1/2 sectors. Making the transformation x — > y x ^ A , ^(x) — > 
y x 4>{y) with A = -(1/2) + 1/{2A), and A = 2K and then rescaling y — > ^foiy with a = {-AK 2 / E) 2K we obtain an 
equation of the same form with parameters if', E given by 

1 1 /AK 2 \ 2K 

<*-*• l' = *-V2, * = —,{—) . (119) 

The point if = 1/2 is a self-dual point of this transformation. Pr esum ably this symmetry is the origin of the 
Seiberg-Witten type duality in the impurity problem observed in Ref. [l39j | . 



K 

(118) 
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2. PT-symmetric quantum mechanics 

The surprising link between the ordinary differential equation (|114[) and the exactly solvable problems related to 
conformal field theories with negative central charge and/or two-color Coulomb gas on a circle is not limited to the 
examples considered above. There is a deep connection between these problems and a non Hermitia n VT -symmetric 
quantum mechanics. The latter has been formulated and studied by C. Bender in various contexts jl4dl |. This story 
goes back to the (unpublished) work of Bessis and Zinn-Justin where, inspired by studies of Lee- Yang model, they 
conjectured that the ID Schrodinger equation with potential V(x) = (ix) 3 has a real spectrum. Later in Ref. |l4l| 
this statement was generalized for arbitrary power-law potentials, V(x) = (ix) n for real n. It was also conjectured 
that the Schrodinger equation is invariant with respect to combined action of the VT symmetry: V (parity inversion): 
x — > — x, p —>■ —p; T : x — ► x, p — > — p, i — ► —i. 

As we have seen, the partition function of the boundary Sine-Gordon model (as well as of the other related models) 
is real for the imaginary value of the coupling constant. This is of course not a generic property of field theory models. 
Therefore the boundary sine-Gordon model with imaginary coupling also belongs to the class of PT-symmetric 
systems. It is then not surprising why the ODE (|1 14[) appears in such a theory. 

This intriguing relation between ordinary differential equations and va rious integrable models has been recently 
extended to a large class of models (for recent review a nd r eferences see |l42j |) by making a close link to the PT- 
symmetric systems. These relationships have been used jl43| | f or th e so-called circular brane model which in certain 
limit describes the so-called Ambegaokar-Eckern-Schon model |l44| . It is therefore reasonable to expect other new 
interesting applications of these findings. 



3. Analysis of distribution functions 

Using exact analytic expressions for Z p (K,ig) at K — 1 after the integral transform we obtain the Poissonian 
distribution for W p (a). The case K — > oo is recovered easily as well. One sees that when p = 0, the distribution 
function is a delta function in the leading order (corrections to this result will be considered in Appendix . As p 
increases, W p (a) rapidly broadens. In the limit of large K, the function W p (a) depends only on the product Kp and 
takes a simple form: 

W p (a) « ( Q / V a > 1. (120) 

This function is peaked at a = 1 for Kp < 1/4, it becomes a step function exactly at Kp =1/4, and for Kp > 1/4. 
W p (a) is a monotonically decreasing function of a. When the product Kp becomes large Kp > 1 the function W p 
becomes Poissonian: W p (a) « 4i£pexp(— AKPa). In general, the tendency of broadening of the distribution function 
remains true for all values of K. 

The distribution given in (|120[) is a particular representative of class of extreme- value distributions called Weibull 
distributions. At p = the limit of K — > oo it will be analytically proven in section IIVD 21 that the normalized 
amplitude a is characterized by the universal Gumbel distribution, for which 

W(a) = Ke K{& - 1) -' 1 - eK{5L ~ 1) ~\ (121) 

where 7 w 0.577 is the Euler gamma-constant. The Gumbel distribution, wh ich also belongs to a class of extreme 
value statistics freq uently appears in various problems (see e.g. boo ks [l24j for his toric introduction), including 
number theory jl45| . 1// noise [lOlj j . Kardar-Parisi-Zhang growth [l46j |. free Bose gas [l47l |. etc. Probably the whole 
distribution function in our interference context can in general be characterized as a certain extreme- value statistics. 



D. Non perturbative solution for the general case 

In this section we present a novel approach for calculating W(a), which is based on a mapping to the statistics of 
random surfaces [HI]. We use this method to evaluate W(a) numerically for a variety of situations. We point out that 
our method is not related to the existence of the exact solution of Sine-Gordon models, but relies only on the structure 
of the multi-point correlation functions in the absence of interactions. As has been discussed in section HVBl W(q;) 
is connected to the partition functions of Sine-Gordon models and the partition functions of Coulomb gases. Our 
mapping provides a new non perturbative tool for calculating partition functions of such models. We expect that there 
should be numerous applications of our new method to other physical problems which can be related to Sine-Gordon 
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and Coulomb gas models. As a concrete application of our approach we prove analytically that for periodic boundary 
conditions in the limit of large K the distri bution of fringe visibilities is given by the Gumbel distribution, one of the 
extreme value statistical distributions Il24ll . 



1. Mapping to the statistics of random surfaces 

We start by observing that G(x,y) is real and symmetric. Hence it can be diagonalized on (0,1) by solving the 
eigenvalue equations 

G(x,y)* f (y)dy = G(m f (x). (122) 

o 

Here / is an integer index, which goes from 1 to oo. "J/ /(x) can be chosen to be real and normalized according to 

l 

* f (x)9 k (x)dx = 5(f,k). (123) 



Then, G(x, y) is given by 

/— oo 

G(x,y)=J2G(f)*f(x)*f(y)- ( 124 ) 
/=i 

Decomposition given by Eqs. ()122j) - (|124() is similar to diagonalization of a symmetric matrix by finding its eigenvectors 
and eigenvalues. Now we can write Zi n from Eq. (|83[) as 



JO 



Z 2n = [ ...J d Ml ...du„d Ul ...^„e^(^'< J G(u -^ )+ ^»< J G(t '-^ ) - 5: M G(u -^ ) ) 

(E; */(«i)) 2 + (Ei */(»<)) 2 -Ei */("i) 2 -Ei */(««) 



JO 



K 



dui...dv n e 

l r i 



-(Ei*/(««))(Ei*/(^)) 



JO 



^...d^e^' ^ */("0-*/(«o) 2 -E:=r(*/(«o 2 +*/(«o 2 )] 



(125) 



Square of the first z— sum in the last line above can be decoupled by introducing Hubbard-Stratonovich integrations 
over auxiliary variables t f : 



l r \ 



/=oc f°° ^ fe -^e E ' t/ V /£ ¥ i (*/( il -)-*/^))-^(*/( M ') 2 +*/( t '-) 2 ) 
Z 2n = I / dux...dv n J] ^7= • (126) 
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Now we can simply integrate over dui, ...,dv n , since all integrals over u— variables are the same (integrals over v- 
variables are also identical): 

z 2n =\n f I <?({*/})"<?({-*/})", (127) 



where 

(128) 



/({*/}) - /"'dx e^/^V^ 1 */^)-^ 1 */^) 2 
Jo 



If all eigenvalues G{f) are negative, then 

&(-{</}) = 9({tf})*, 9({tf})g({-t f }) = \g({t f })\ 2 . 
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FIG. 10: Scaled distribution functions of the normalized interference amplitude W(a) at T — for ID gases with open boundary 
conditions, shown for Luttinger parameters K = 2 (dashed), K = 3 (dotted), and K = 5 (solid). The inset shows a comparison 
between open (solid) and periodic (dashed) boundary conditions for K — 5. The figure is taken from Ref. [581 ] . 



^From comparison of Eqs. ([55]) and (|127p we obtain the central result of this section 

t 2 

/=°° r°° e —£dtf 

W(a) = J] J "°° ^ [a - | fl ({t,})| 2 ] . (129) 

/= i V^vr 

Eq. (|129p can be used to simulate distributions W(a) using Monte-Carlo approach. First, one needs to solve integral 
Eqs. (|122p - p24|) numerically to obtain eigen functions and eigen vectors. Then one needs to choose random numbers 
{tf} from the Gaussian ensemble, and plot the histogram of the results for |<?({t/})| 2 . In what follows we perform 
simulations of W(a) with up to N = 10 6 — 10 7 realizations of {tf} and smooth the data. We use a finite value of f max 
and check for convergence with f ma x, typically ~ 30. (a) is always kept within 1% from its expected value. For most 
of the presented results, all eigenvalues G(f) are negative, and Eq. (|129[) can be directly applied. Special care should 
be taken for ID case with nonzero temperature, where one of the eigenvalues G(l) can be positive. This situation 
can be handled by subtracting sufficiently large positive constant C from G(x,y, which makes all eigenvalues 
negative. According to (|53"]) . this leads to rescaling of a by a factor e~ c , which can be easily taken into account. 

In Fig. [TO] we show scaled distribution of contrast VF(<5) at T = for ID gases with OBC for various K. The inset 
shows a comparison between OBC and PBC for K = 5. In Fig. [11] we show the scaled distribution of contrast for 
ID gas with OBC at nonzero temperature and K = 5. As has been discussed earlier, for ^tK/L <C 1 distribution is 
Poissonian and wide, while for K ^S> 1 and £tK/L ^ 1 it is very narrow. Evolution of the full distribution function of 
the visibilities while L is varied can be used to measure the thermal length , £y — hv s / '(fcsT), precisely and to extract 
the temperature. As seen in Fig. [TT] at T ^ the distribution function has characteristic features, i.e. it is generally 
non-symmetric and can have a minimum. These features can be used to distinguish the noise due to fluctuations of 
the phase from technical noise. Finally, in Fig. [TO] we show scaled distribution of contrast for 2D gas with unity aspect 
ratio of imaging area and OBC below BKT temperature. Above BKT temperature distribution functions become 
Poissonian for L ^> £, where £ is the correlation length. In 2D one cannot describe the crossover at L ~ £ similar to 
ID, since the action which describes the fluctuations of the phase is not quadratic in this region, and Eq. I|83p doesn't 
hold. 

Result (|129|) can be interpreted as the mapping of our problem to the statistics of random surfaces (strings in ID), 
subject to classical noise. In this mapping, *B f(x) are the eigen modes of the surface vibrations, tf are the fluctuating 
mode amplitudes, and |G(/)| is the noise power. This mapping holds as long as all eigenvalues G(f) are negative, as 
discussed earlier. Infinite dimensional integral over {tf} variables can be understood as an averaging over fluctuations 
of the surface. For particular realization of noise variables {tf}, complex valued surface coordinate at point x is given 
by 



/ 



Kx; {t f }) = J2 tf\/^P-*f(z) ~ ^*/(*) a - ( 13 °) 



g({tf}) is a number defined for each realization of a random surface {tf}, given by Eq. (|128|) . Fringe amplitude a for 
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FIG. 11: Scaled distribution functions of the normalized interference amplitude W(q) for a ID Bose gas with open boundary 
conditions at nonzero temperature and K = 5. Different curves correspond to ratios K^t/L = oo (solid), K^t/L = 1 (dotted), 
and K^t/L = 0.25 (dashed), fr is the thermal correlation length, K is the Luttinger parameter and L is the imaging length. 
The figure is taken from Ref. [581 ] . 




FIG. 12: Scaled distribution functions of the normalized interference amplitude W(a) for a two dimensional Bose gas with 
the aspect ratio of the imaging area equal to unity and open boundary conditions. Temperature is below the Berezinskii- 
Kosterlitz-Thouless (BKT) transition temperature. Different curves correspond to T = Tbkt ( BKT transition point, solid), 
T — (2/3)Tbkt (dashed line) and T = (2/5)Tbkt (dotted line). Above the BKT transition temperature the distribution 
function is Poissonian. The figure is taken from Ref. [5q . 



each realization of {£/} variables is given by 



= \9({tf})\ 2 = 



dxe H^{t f }) 



(131) 



2. From interference of ID Bose liquids of weakly interacting atoms to extreme value statistics 

To illustrate the power of the interpretation of the interference fringe amplitudes in terms of random surfaces, 
we now prove analytically that for periodic boundary conditions in ID and in the limit of large K, the normalized 
distribution W(a) is given by the Gumbel distribution (|12ip . one of the extreme value statistical distributions. We 
note that to prove such result using 1/K expansion described in Appendix [5] one needs to go to the infinite order of 
the perturbation theory, so essentially our result is non-perturbative in K. 

The Gumbel function 

P G (X) = e -(*+7)-e-<*+^ (132) 



plays the same role in the extreme value statistics 124| as the usual Gaussian distribution plays in the statistics of 
the average value. According to the central limit theorem, the average value of numbers taken from the same 
ensemble in the limit of large Af is distributed according to Gaussian function. One can prove similar theorems for 
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the distribution of the extreme, i.e. largest value of N 1 numbers taken from the same ensemble. Gumbel function 
(|132[) is one of the universal functions describing ensembles which decay faster than any algebraic function at infinity. 
The Gumbel distribution often appears in applied mathematics, including problems in finance or climate studies, 
where it is used to describe rare events such as stock market crashes or earthquakes. 

For the periodic case corresponding to Eq. (|85|) , the eigenvalue problem (|122|) - (jl24|) can be solved analytically, and 
the noise has 1// power spectrum: 

1 /= °° 1 r 

G per (x,y) =log-sm7r|x-2/| = -log27r- V — (y/2 cos 2irfx)(V2 cos 2-jrfy) + (V2mx2nfx) (y/2 sin 2irfy) 

7T ^ — ' 2 f L 

/=1 J 

(133) 

Eigen modes are given by simple harmonic functions, and all of them, except one, are doubly degenerate: we will 
denote corresponding noise variables as t\j,t2j and eigen values as G(f) = —1/(2/). Simplification in the limit 
K ^> 1 stems from the fact that exponent in Eq. (|131[) can be expanded in Taylor series, since h(x; {tf}) is small for 
large K according to its definition Q130P . In this case it can be shown th at th e distribution of a is linearly relat ed to 
roughness, or mean square fluctuation of the surface, as defined in Ref. [lOl| . It has been shown in Ref. [lOl| that 
1/ f noise results in th e Gumbel distribution of the roughness, which has been interpreted in terms of extreme value 
statistics in Ref. [l48j |. 

Terms of the order of 1/yK vanish in the Taylor expansion of Eq. (|13ip . since the average values of cos (2tt fx) 
and sin(27r/a;) on interval (0, 1) are equal to 0. To order 1/K, we obtain 

a « 1 + 1 10S27T - f) JL(t? (/ + tl <} - 2). (134) 
/=i J 

The constant term — log27r in G per (x, y) gives a constant rescaling of the distribution function of a, and doesn't show 
up in the normalized fringe amplitude a : 

oo ^ 

5 « 1 -E 2 T^,/+4/- 2 )- (135) 
Thus to the leading order in 1/K expansion, the distribution Y(x) of the rescaled variable 

x = -K(a-l) = jr±;(tl f +4 tf -2) (136) 



equals 



/=0O „oo „oo , , 2 ,,2 / /=00 



f=l " » -«j v /=1 



To prove Eq. (|121[) . we need to show Y(x) — Pc(x), where Gumbel function Pg(x) is given by Eq. (|132|) . Let us 
introduce new positive variables 

«/ = J 2/ J > 0, (138) 

then 

/=°o /.OO . / /=°° 



J poo I ~ 



(139) 



To prove that y (a;) = Pg{ x ) w e will calculate their Fourier transforms: 



dxY{x)e lsx = Yl fdu f e- fu ie ls[u t- 1/f) = f[ 1 ' . = e - 47S r[l - is], (140) 

"°° /=i /=i ^ * S 

/•OO /»OC 

P G {is)= I dxP G {x)e lsx = dxe- {x+ ^- e e isx = e-^ s T[l-is}. (141) 
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The proof abov e doe sn't illustrate the meaning of Y(x) as of a distribution of extreme value. Here we will follow the 
method of Ref. jl48| | and explicitly construct the variable, extreme value of which generates Gumbel distribution. Let 
us impose a finite cutoff f m ax = N, and at the end of calculation we will send N to infinity. If one identifies 

z\ = u N (142) 

z 2 = un-i + u N , (143) 

(144) 

zn = u\ + u 2 + ... + u N , (145) 
then {zi, z^} is an ordered set (since Ui > 0) of outcomes taken from Poissonian distribution, since 

T f=1 /»/ _ e -EJ=i *f } (146) 
and Jacobian of transformation from variables {«/} to {zf} variables is unity. Then 

Y N (x)=N\ f e~ Zl d Zl [ e- Z2 dz 2 ... [ e~ ZN dz N 5 [ x - ( z N - ^ 4 ) J (147) 

J0 J z\ Jzn-i \ \ t—i J J J 

is nothing but the shifted distribution of the largest of N numbers taken from the Poissonian distribution, and in the 
limit of large N this distribution converges to Gumbel function. 

One can understand the appearance of the Gumbel distribution by noting that for K 3> 1 the distribution function 
of the interference amplitude is dominated by rare fluctuations of the random periodic ID string, which are spatially 
well localized. The Gumbel distribution was introduced precisely to describe similar rare events such as stock market 
crashes or earthquakes. For open boundary conditions the u niversal distribution for large K is slightly different from 
Gumbel function, similar to 1// noise in other systems |l0l| . But the main properties, like the presence of asymmetry 
or the asymptotic form of the tails are preserved. 



V. CONCLUSIONS 



A. Summary 

When we discuss interference experiments with ultracold atoms, the conventional idea of the particle-wave duality 
takes a new meaning. On the one hand, these experiments probe phase coherence which is typically associated with 
coherent non-interacting waves. On the other hand, one can use powerful tools of atomic physics to control interactions 
between atoms in a wide range and to reach the regime of strong correlations. One can also prepare atomic systems 
in states which would be difficult if not impossible to obtain in optics, e.g. low dimensional condensates with strong 
thermal or quantum fluctuations. This remarkable combination places interference experiments with ultracold atoms 
in a unique position: they can address a problem of how the interactions, correlations, and fluctuations affect the 
coherent properties of matter. This question appears in many areas of physics, including high energy and condensed 
matter physics, nonlinear quantum optics, and quantum information. While the naive answer that interactions 
suppress interference turns out to be correct in most cases, the goal of these lecture notes was to demonstrate that 
the quantitative analysis of this suppression can provide a lot of nontrivial information about the original systems. 

We discussed two effects which contribute to the reduction of the interference fringe contrast in matter interfer- 
ometers. The first effect is the shot noise arising from a finite number of atoms used in a single measurement. This 
analysis is particularly important for interference experiments with independent condensates in which the position 
of interference fringes is random and averaging over many shots can not be performed. In this case one needs to 
rely on single shot measurements to observe interference patterns. While interference of independent condensates has 
been discussed before [6(1 [6l|, [62|, HH, [H, IsS 0, [5E [sE IIzL Isl], to our knowledge, we provide the first derivation 
of the full distribution function of the amplitude of interference fringes. Another mechanism of the suppression of 
the amplitude of interference fringes discussed in these lecture notes is the quantum and thermal fluctuations of the 
order parameter in low dimensional condensates. The motivation for this discussion comes from the observation that 
interference experiments between independent fluctuating condensates can be used to study correlation functions in 
such systems [55] ■ For example, one can use the scaling of the integrated amplitude of interference patterns to analyze 
two point correlation functions. This method has been successfully applied by Hadzibabic et al. [12] to observe the 
Berezinskii-Kosterlitz-Thouless transition in two dimensional condensates. One conceptual approach to understanding 
interference experiments with independent condensates is to consider them as analogues of the Hanbury Brown and 
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Twiss experiments in optics [7]. In the latter experiments interference between incoherent light sources appears not 
in the average signal but in the higher order correlation function. One important difference however is that matter 
interference experiments are of a single shot type and information is contained not only in the average fringe contrast 
but also in the variation of the signal between individual shots. In particular higher moments of the amplitude of inter- 
ference fringes contain information about higher order correlation functions [5a ] . A complete theoretical description of 
the fringe contrast variations is contained in the full distribution functions of the fringe amplitudes, which we calculate 
for one and two dimensional condensates [56l . |58| in the limit when the number of atoms is large and the shot noise 
can be neglected [99j . An important aspect of these lecture notes was identifying intriguing mathematical connections 
which exist between the problem of calculating distribution functions of interference fringe a mpli tudes and several 
other problems in field the ory a nd statistical physics, such as the quantum impurit y problem [100(, tunneling in the 
presence of the dissipation [129], non-hermitian PT-symmetric quantum mechanics [l40l Il4ij and various conformal 
field theories. We developed a novel mapping of a wide class of such problems to the statistics of random surfaces, 
which provided a complete non perturbative solution. In certain cases we have analytically proven (58j the relati onsh ip 
between the distribution function of fringe amplitudes and the universal extreme value statistical distribution [l24j . 

B. Some experimental issues 

We now comment on a few issues relevant for experimental analysis of noise in interference experiments. The amount 
of information contained in the experimentally measured distribution function is directly related to the number of 
cumulants which can be accurately extracted. This includes the second cumulant ki , which corresponds to the width 
of the distribution; the third cumulant &3, which is related to skewness, g\ — k^/k^ 2 , and describes the asymmetry 
of the distribution function, and so on. In general, the statistical error in determining the n th order cumulant after 
N measurements scales as \J A n /N, where A n is a constant which grows with n and depends on the higher moments 
of the distribution. For example, to experimentally distinguish the normal and Gumbel distributions it is necessary 
that the statistical error in skewness is at least a factor of two smaller than the mean skewness, which for the Gumbel 
distribution is g\ » —1.14. Thus the minimal number of measurements required is N m i n sa 24/^ w 20, where we 
used A% sa 6, appropriate for the normal distribution [l49j |. In practice the required number of measurements may be 
higher because of the influence of other possible sources of noise. However, it is certainly experimentally feasible. 

Another experimentally relevant issue is the effect of the inhomogeneous density due to the parabolic confining 
potential. While the approach discussed in these lecture notes can be extended to include the inhomogeneous density 
profile, interpretation of the experimental results is more straightforward when density variations can be neglected. 
We note that when the c onde nsate density varies gradually in space, the power-law decay of the correlation functions 
is not strongly affected jl50l ]. except that the exponent may be different in different parts of the trap (correlation 
function exponents typically depend on the density) . We expect that qualitatively our results will not change provided 
that the power law decay of the correlation functions is much stronger than the change of the condensate density in the 
measured part of the cloud. To be more precise, the best comparison with theory can be achieved when the obser vation 
region L is much smaller than the size of the condensate, determined by the effective Thomas- Fermi length [l5lj j . 
Rtf — {SNT1 2 1 ' (m 2 uj 2 aiD)) 1 ^ (here N is the number of atoms of mass m and am is the one-dimensional scattering 
length). As long as L remains much larger the healing length, our analysis is valid. In the regime of weakly interacting 
atoms, one can show that the ratio between the effective Luttinger parameter K at the center of the harmonic trap 
and at the boundary of the observation region is given by 1 — L 2 /8R^ F , thus giving only a small correction to the 
distribution function computed in the central region. One can also reach similar conclusions in the strongly interacting 
regime. It is also worth pointing out that we expect the limiting case of the Poissonian distribution to be particularly 
robust to the inhomogeneous density of atoms. Indeed the Poissonian distribution is related to the fast 1/y/x decay 
of the one-particle correlation functions in the strongly interacting limit. This scaling is a universal feature of the 
Tonks-Girardeau limit of bosons and is not affected by the weak harmonic trap [152]. 

C. Outlook 

Before concluding these lecture notes we would like to discuss questions which still need to be understood in the 
context of interference experiments with ultracold atoms. We also suggest an outlook for future theoretical work. 

Combining shot noise with the order parameter fluctuations. 

A careful reader has undoubtedly noticed that we discussed either the shot noise or the order parameter fluctuations. 
At this point we are still lacking theoretical tools which would allow to include both effects simultaneously. One of 
the difficulties is that such analysis requires the knowledge of the correlation functions for all distances rather than 
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the long distance asymptotic form. Indeed, in section IIVI we showed' that the short distance part of the correlation 
functions gives contribution of the same order as the shot noise. In the particular case of the interference of 2D 
condensates, the knowledge of the short distance behavior of the correlation functions is needed to include the effect 
of the vortex excitations below the BKT transition. 

Stacks of independent condensates 

In these lecture notes we focused on interference patterns from a single pair of condensates. However in experiments 
one often has a stack of several condensates (see e.g. [3l|). In this case interference arises from all possible pairs, and 
the system provides intrinsic averaging and suppression of the noise. For a finite number of condensates self-averaging 
is not complete and one expects finite fluctuations of the fringe contrast. It would be useful to generalize analysis of 
the shot noise and order parameter fluctuations to such systems. 

Dynamics of interacting atoms. 

One of the advantages of the cold atoms systems is the possibility to study non equilibrium coherent dynamics of 
interacting systems. In particular dynamical splitting of a single condensate into a pair of condensates has been 
performed in expe riments on microchips [13, EH, |H, |4^| and stimulated theoretical work on the subject jl53l Il54[ 
Il55l . 1 156L Il57l Il58| . Similar e xper i men ts can also be done using superlattice potentials in optical lattices which are 
now available in experiments jl59l Il60l ] . While analysis of fringe amplitude distribution functions presented in these 
lecture notes dealt exclusively with systems in the thermodynamic equilibrium, it would be interesting to generalize 
it to systems undergoing nonequilibrium dynamical evolution. 

Interference experiments with fermions 

The discussion present ed in these notes was limited to the case of interference of boson s. Such exper iments can also 
be done with fermions [l6l| . which are available experimentally in diffe rent dime nsions [162L Il63l Il64| . For fermions, 
modulation of the density can be related to fermion antibunching 20, 165, 166]. Analysis of the noise of the fringe 
contrast visibility for fermions would be an interesting problem too. 

Generalization to other systems 

We note that mapping of the Coulomb gas into the statistics of random surfaces introduced in section IIVDI should 
have applications beyond calculating the distribution functions of the interference fringe amplitudes. This is a new non 
perturbative tool to calculate partition functions of a variety of other systems that can be represented as Coulomb gas 
models. Examples include quantum impurity-related problems [l00| , Sine-Gordon models where interaction strength 
can depend on position, and many others. Our mapping is not related to the existence of the exact solution of Sine- 
Gordon models, but relies only on the factorable structure of the many-point correlation functions in the absence of 
interactions, which is a general property of a Gaussian action. 
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APPENDIX A: LARGE K EXPANSION 

In this Appendix, we will describe a systematic " diagrammatic technique" to calculate Zi n or as an expansion 
in small parameter 1/K. It corresponds to the "high temperature" limit of the classical gas analogy discussed in 
section ITV Al This expansion can be applied both in ID or 2D, and can be used to study the limiting di strib ution at 
large K, which for PBC in ID has been conjectured [5(| and proven [58| to be the Gumbel distribution [l24j . 

1. Expansion to order (1/K) 2 . 
We will start from the ID case by expanding the exponent in Eq. 1831) : 
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G(x,y) G(x,y) G ( X; y) 

x y y x x y 

G(z,t) I G (y,z) | G(x,y) | 

z t y z x y 



FIG. 13: Topologically inequivalent diagrams, corresponding to 1/K 2 terms in expansion of Zi n 

Z 2n = / ... / du 1 ..Au n dvi..Av n e^^<o G{u ^^ )+ ^ G{v - v ^-^ G[u ^^ ) ) = 
Jo Jo 

I ■■■ I du 1 ...du n dv 1 ...dv n ^ — — I G(u t , Uj) + ^ G(vj, vj) - G(u i: vj) 

J° J° m=0 m ' \i<j i<j ij j 

In the first order of expansion in powers of 1/K, G(x,y) dependence comes only through one integral 



(Al) 



1 r l 



I Q = / dxdyG(x,y) (A2) 
Jo Jo 

after the integration in (|A1[) . 

The prefactor depends only on n, and can be calculated analytically: the total number of u — u terms is n(n— l)/2, 
the total number of v — v terms is also n(n — l)/2, and the total number of u — v terms is n 2 . The latter come with 
— 1 sign, so the total expression for Z 2n up to 0(l/K 2 ) is 

Z 2n = l-±nI +O( 1 ^). (A3) 
In second order, dependence of Z 2n on G{x, y) comes through three different integrals: 

l r i r i r i / r i r i \ 2 

2 



h,2= / / / dxdydzdtG{x,y)G{z,t) = / / dxdyG{x,y)\ = l£, (A4) 
Jo Jo Jo Jo \Jo Jo / 

h,2= / / / dxdydzG(x,y)G(y,z), (A5) 
Jo Jo Jo 

h,2= I I dxdyG(x,y)G(x,y). (A6) 
Jo Jo 

Pictorially, these expressions can be represented by the corresponding diagrams, shown in Fig. 1131 There, horizontal 
solid line corresponds to G(xi,Xj). If two ends are connected by a dashed line, then in the integral these two ends 
should correspond to the same variable. To write the expression for Z 2n in 1/K 2 order, one has to calculate the total 
number of terms corresponding to Ii >2 , 1 2 . 2 , I^^ 2 . Clearly, it will be some universal polynomial depending on n, which 
is determined by combinatorics. Here we describe how to do this combinatorics in detail. 

When parenthesis are multiplied in (|Alj) there are several types of expressions which appear, and come with different 
signs to the integral: 



ua,Vji)G(u i2 ,Vj 2 ) + 

E (G(uu,Uji)G(u i2 ,Uj 2 ) + G(v a ,Vji)G(v i2 ,Vj 2 ) + 2G(u ll ,u jl )G(v l2 ,v j2 )) 

il<jl,i 2 <j 2 

-2 ^ (G ! ('Uii,%i)G ! (u i 2,«i2) + G(va,v jl )G(u i2 ,v j2 )). (A7) 

il<jl,i 2 ,j 2 



3G 



A) B) 

FIG. 14: Topologically equivalent diagrams, corresponding to a)/2,2 and b)/3,2. 



These are (signs are indicated) 



+ G(uu, Uji)G(ui2, Uj2), and equivalent after integration G(vn, Vji)G(vi2, 1^2); 

+G(Uil,Vjl)G(Ui2, Vj2 ) j 

-G(un, Uji)G(v,i2, Vj'2), and equivalent after integration — G(vn, Vji)G(ui2, Vj2); 

+G(uii,Uji)G(vi2,Vj2); 



(A8) 
(A9) 
(A10) 
(All) 



For terms which have G(v,ik, Ujk) or G{vik 1 Vjk), k = {1, 2} there is a restriction that ik < jk, but there is no such 
restriction for G{uik, Vjk)- 

Lets calculate the total number of 2i l2 , h,2 and ^3,2 terms for u — u — u— u expressions of type (|A8j) . Total number 
of /i : 2 terms compatible with restrictions is 



^n(n - l)(n-2)(n-3). 



(A12) 



1/2! 2 appears since pairs and {i2, j2} are ordered , while n(n — l)(n — 2)(n — 3) is the total number of ways 

to choose a sequence of four different not ordered numbers out of the set of size n. 
Total number of ^2,2 terms compatible with restrictions is 



2! 2 



n(n - l)(n - 2) * 4. 



(A13) 



Here the situation is similar to I12, but there is a "topological" prefactor of 4, which corresponds to four topologically 
equivalent diagrams for 7 2 ,2, as shown in Fig. [T4la. 
Finally, total number of 7^2 terms is 



1 

2!2 



n(n — 1) * 2. 



(A14) 



Here, 2 is again a "topological" prefactor, corresponding to two topologically equivalent diagrams shown in Fig. fT4lb. 
Overall, u — u — u — u term of (IA8D gives 



-n(n - l)(n - 2)(n - 3)J li2 + n(n - l)(n - 2)/ 2 , 2 + -n(n - 1)Z 3 , 2 (A15) 
As a simple check of combinatorics one can calculate the total numbers otu — u — u — u terms of type (|A8|) . It is 



1 n(n- l)(n- 2)(n-3) + ^n(n - l)(n - 2) * 4 + ^n(n - 1) * 2 = " "'" ' 



2! 2 



2! 2 



2! 



(A16) 



as it should be from the calculation which neglects the dashed lines, and doesn't impose any restrictions on the terms 
at different horizontal rows. 

Analogously, one can calculate polynomials for expressions (IA9|) - (IA11|) . and results are as follows, respectively: 



n 2 (n - l) 2 / li2 + 2n 2 (n - 1)J 2;2 + n 2 / 3i2 , 
-\n 2 {n - l)(n - 2)I h2 - n 2 (n - l)/ 2 , 2 , 



2 2 



n 2 {n-lfl x< 



(A17) 
(A18) 

(A19) 
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Summing all the terms with corresponding prefactors from (|A7[) . we obtain 

Z2n = 1 - -j-nlo + — 2 (Pl )2 /l,2 + P 2 ,2h,2 + P^hfi) + °(^) = 

1 - ^rd + ^ (3n(n - l)J li2 - 4n(n - l)/ 2 , 2 + (2n - l)nJ 3 , 2 ) + 0(^3), (A2Q) 

where Pi :2 — 3n(n — 1), P 2 . 2 = — 4n(n — 1), and P3.2 = (2n — l)n are universal polynomials of n. Notice, that due to 
sign cancellations of terms in (|A8|) — (I Al 1|) . the overall degree of polynomials in 1/K 2 order is m = 2, compared to 
naively expected 2m = 4. 

Using (|A20p . one can calculate Z 2n /Z% up to (9(l/iv 3 ) : 

if = 1 + !! ^ ! ^( / ^ - 27 2i2 + 7 3)2 ) + 0(^3) (A21) 

Calculated values of {J , h,2,h,2, h, 2 } an( A {PcT^ ^1 T' ^2 jm ^32 } are 

3 



7 = -| = -1.5, 4 ,er = -log27r w -1.83788 
Ji,2 = \ = 2.25, JfJ = (log27r) 2 « 3.37779 

7 2 .2 = 51 7/ 2 » 2.28502, If 2 r = (log27r) 2 « 3.37779 
18 



(A22) 
(A23) 
(A24) 



= J - 3.5, Iff - - J + y' « 4.20026. (A25) 



Thus, for if — > 00, the limiting ratio of the widths of the distributions for OBC and PBC are 



Z A /(Z 2 ) 2 -1 h,i-2h, 2 + l. 



3,2 



zr/(z p 2 er ) 2 - 1 (i p j - 2/2T + 11%) 



1.43465. (A26) 



2. General properties of (1/K) m terms, and expansion to order (l/K) 5 . 

^From the expansion of the previous subsection, one can formulate the general properties of the "diagram technique" 
to calculate terms up to (l/K) m : 

a) First, one has to draw all possible topologically inequivalent diagrams, which consist of m horizontal solid lines, 
with some of the ends connected by dashed lines. Each end can have at most two dashed lines coming out of it. 

b) Ends which are connected by a dashed line correspond to the same variable. Diagrams for which two opposite 
ends of the horizontal line correspond to the same variable are excluded. 

c) Expression which corresponds to a diagram is constructed the following way: if variables at the end of a given 
solid horizontal line are x and y, then G(x,y) should be put as one of the terms in the product under the 
integrand. Thus the integrand consists of the product of function G of some variables m times. Diagrams for 
which the integrands are the same up to relabeling of the variables are considered to be identical. 

d) all free variables should be integrated from to 1. 

Diagrams can be connected or disconnected. For example, in Fig. 1 131 diagram corresponding to I1.2 is disconnected, 
and diagrams corresponding to Z 2>2 and 1^^ are connected. The integral which corresponds to a disconnected diagram 
is a product of expressions, corresponding to its parts: for example, Ji )2 = Io* Io- 

If the number of topologically inequivalent diagrams of order m is g(m), then the term of the order 1/K m in Z 2n 
has the following form: 



r 1 1 ( \ " 1 9{m) 

. d Ul ...dv n ^ Km ^G(v,i,Uj) + ^2G(vi,Vj) - ^2G(ui,Vj) = m]Km (^2 Pr,m(n)J r , m ), (A27) 
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where P r ,m(n) are universal polynomials of n, which can be calculated combinatorially, as described in section [A 11 
Polynomials P r<m (n) should satisfy the following requirements: 

a) For positive integer n, the values of P r ,m(n) are integer, and 

P r>m (0) = 0. (A28) 



b) If one sets G(x, y) = 1, then I r m — 1, and lhs of (|A27|) can be trivially calculated. This implies 

aim.) 

p -»» = (- n ) m - ( A29 ) 

r=l 



c) Degree of each polynomial P r ,m{n) is not larger than m. This has been shown above for m — 1 and m — 2, and 
has been checked up to m = 5, although we didn't succeed in proving it directly. This conjecture is supported 
by the fact, that it guarantees that for any G(x,y) for K — > oo distributions fit on the top of each other after 
proper rescaling (see next section). 

One can do the combinatorial calculations similar to previous section for m — 3, and it takes about a day "by 
hand". We will only show the results here. For m — 3 there are ^(3) = 8 topologically inequivalcnt diagrams, which 
are shown in Fig. 1151 out of which 5 diagrams (Ii,3 to 15,3) are connected. Expressions corresponding to each diagram 
and universal polynomials are, respectively: 



f.1 r l fl el 

h,. 



JO JO JO 

W/77' 

Jo Jo Jo Jo 

fl /-l /-l 

13,3 = 



Jo Jo 



'4,3 



1 rl r l 



Jo Jo 



'5,3 



- l)(2n- 


3), 


(A30) 


12n(n — 


1), 


(A31) 


-12n(n- 


1), 


(A32) 


- l)(2n- 


1), 


(A33) 


5,3 (n) = - 


-n, 


(A34) 


-l)(n- 


2), 


(A35) 


- l)(2n- 


3), 


(A36) 


-l)(n- 


2). 


(A37) 



Numerically evaluated integrals for OBC and PBC are 



7 13 = -3.49399, If J = -6.20797; (A38) 

h,3 = -3.5268, JfJ = -6.20797; (A39) 

Z 3)3 = -5.32704, ZfJ = -7.71956; (A40) 

h, 3 = -4.21255, 7|J = -6.50848; (A41) 

I 5i3 = -11.25, 7fJ = -12.5458; (A42) 

h,3 = -3.42753, TfJ = -6.20797; (A43) 

J 7 ,3 = -5.25, J?J = -7.71956; (A44) 

7 8)3 = -3.375, ifj = -6.20797. (A45) 

For m > 3 it becomes too cumbersome to manually calculate universal polynomials P r ,m- We wrote a program in 



Mathematical which expands m — th term of (|A1|) directly, and calculates the values {P r>m (0), P r ,m( n )} using 
powerful pattern recognition tools. After that, P,-. m (n) is recovered using Newton's formula. Results for m — 3 can 
be recovered that way. One can also check in each order that the degree of P r ,m{n) is not larger than m. For each m 
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FIG. 15: Topologically inequivalent diagrams for m — 3. Diagrams Ii t s — Is,s are connected and — Ia,3 are disconnected. 



FIG. 16: Irreducible diagrams for m = 4. 



the program needs as an input all topologically inequivalent diagrams, and currently the results have been extended 
to m = 5. For m = 4, the overall number of diagrams is 23, out of which 12 are irreducible, and shown in Fig. 1161 
For to = 5, the overall number of diagrams is 66, out of which 33 are irreducible. Numerical prefactors of polynomials 
Pr.m(n) grow with m, while their overall sum has a prefactor 1, as follows from (|A29[) . For example, one of the 
polynomials for m = 5 has a prefactor of 384. This puts stringent requirements on the errors in calculation of I r ^ m {n), 
so going beyond m = 5 will require additional numerical effort. For to = 5 results are reliable, since we can check the 
numerical accuracy of the calculation of integrals by comparison with analytically proven [5a ] Gumbel distribution. 

While calculation of original Zi n requires In — dimensional numerical integration, calculation of the I r ^ m requires 
at most to — dimensional numerical integration. This can be seen from the fact that even for irreducible diagrams, 
only the parts which contain "loops" have to be integrated numerically. Horizontal lines which have a free end can 
be "integrated out" analytically, and the dimension of the numerical integral has at most dimension of the loop. For 
example, first and third horizontal bars in ii 3 of Fig [15] can be integrated out using the following identity: 



dy log \x — y\ = — 1 + log (1 — x) — a; log (1 — x) + xlogx, 



(A46) 



so one has to do only 2-dimensional integral numerically. Analogously, diagram ^2,3 requires only one dimensional 
numerical integration. In each order m there is only one diagram which requires m — dimensional integration, i.e. I4.3 
in Fig |15l for m = 3. All the rest require at most (to — 1) — dimensional integration. 
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3. Properties of the K — > oo distribution 



For K — ► oo distribution function becomes very narrow, and it is an interesting question to investigate the limiting 
behavior of the distribution function. Let us consider the distribution of normalized contrast a, defined by (|88p . Due 
to normalization the following relations hold: 



W{a)da = 1, / aW(a)d& = 1. 
ii Jo 



(A47) 



For large K distribution function is peaked near 5 = 1, and its width is proportional to 1/K. To calculate the 
properties of the universal function we assume that it decays relatively fast (exponential decay is enough) away from 
a = l for large K, so we can extend integrations to infinity. If we define the fluctuation of the normalized contrast 
(3 = 5 — 1 with distribution function Wo{p) = W{1 + (3), then 



- 1 



£2 
Z 2 



{a - l)W{a)da 



11 



-1 



Z, 



z & 1 , 



-1 + -J = / {of- l)W(&)da 

poo 



l)W(5)da 



(2/3 + /3 2 )W (/3)d/3 = 2Afi + M a , 
(3/3 + 3/3^ + /3 3 )WV/?)d/3 = 3Mi + 3M 2 + M 3 , 



(A48) 

(A49) 

(A50) 
(A51) 



where Mi, M2, M3 are the moments of Wq{$). Inverting (|A48|) - (|A50|) . one can find M2, M3, ... as series expansion in 
powers of 1/K. If for large K the distribution functions collapse on the top of each other after proper rescaling, as 
has been conjectured in [5(|, then the expansion of M m in powers of 1/K should start from 1/K m . Below we will 
show, that this is necessarily true for all n, if the degrees of universal polynomials P r m (n) are not larger than m and 
P r , m (0) = 0. 

To show this, let us consider an analog of Eqs. (|A48|) - (|A50[) for the fluctuation of the unnormalized contrast 
P = a — 1 with distribution Wq{(3) — W{1 + (3) (normalized contrast considered in Appendix IA 31 is related to a by 
simple rescaling). Then 



a n W{a)da 



{l + f3) n W {(3)dp 



(A52) 



i=0 



i=0 



where Mj is the i — th moment of Wq{0). We need to show that expansion of Mj in powers of 1/K starts only from 
the terms of the order of 1/K 1 . This can be seen from expansion (|A27[) together with (|A52[) . If degrees of P r ,m(«) are 
not larger than to, then (|A27[) means that 1/K m terms in Zm grow at most as n m for large n. On the other hand, 
C« — Tt "i y grows as n l for large n, and if expansion of Mi in powers of 1/K starts before 1/K\ this will contradict 



i\(n— i)\ 

the previous sentence. 

To find the first nontrivial contribution to M TO , one has to go to order 1/K m in the expansion. Using results 
for tyi — 5 calculated above, one can calculate for K — > 00 limiting behavior of K 2 M2, ...,K 5 M$ for periodic and 
non-periodic cases, and compare it with the result of the Gumbel distribution 



Wq0) = Ke K ^- e 
where 7 w 0.577 is the Euler gamma-constant. One obtains: 



Non-periodic : K 2 M 2 2.35991, K 3 M 3 



-5.105577, K^JVU -> 37.5258, K 5 M 5 



Periodic : K 2 M% e 



1.64493, K 3 Mg er -> -2.40411, K 4 M% e 



Gumbel : K 2 M? 



1.64493, K 3 M^ 



-2.40411, K 4 M? 



14.6114, K 5 MP 



M 2 

-242.492, -J- 
M| 

14.6114, K b Ml er — s 

^_ 64 .4321,M^ 

(Mf )3 



(A53) 

1.98336(A54) 
-64.4321(A55) 
> 1.29857(A56) 
(A57) 
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Clearly, for periodic case the agreement with the Gumbel distribution is excellent, but for non-periodic case, one can 
unambiguously conclude that limiting function is NOT a Gumbel function. For non-periodic case the function is more 
widely distributed compared to the Gumbel function, as evident from the higher moments. Here for completeness we 
provide numerical results for Zzn , Z^ up to m = 5 : 

1.8379/1 -0.4112n + 2.5114?i 2 0.1002n - 0.1548n 2 + 2.1456n 3 
Z 2n - H I ^2 1 ^3 l ~ 

-0.03382n + 1.0804n 2 - 0.7399n 3 + 1.7369/x 4 0.01535n + 1.05941ti 2 - 0.161727n 3 + 0.172302n 4 + 0.93098n 5 



K A K 5 

1.5n -0.554956n + 2.30496n 2 0.049657n + 0.343838n 2 + 1.481505n 3 
Z 2n = 1 + — + W2 + ^ + 

-0.04741n + 2.0329n 2 - 1.8735?i 3 + 1.8256n 4 0.0106468n + 3.14868n 2 - 4.1829n 3 + 2.8980n 4 + 0.0943039n 5 
K 1 + K~ 5 ■ 



4. D=2 



Here we will briefly consider results of the expansion in powers of 1/K for 2-dimensional case. Correlation function 
below the BKT transition is given by Eq. (|78[) . and for square imaging area with unity aspect ratio, all discussions of 
one dimensional case carry over, with substitutions 



J rl i-l pi pi pi 

2„-? / J„. . / / j2 



dui -> / / d Ui, dvi^ / d Vi,G(x,y) = log \x - y\, (A58) 
Jo Jo Jo Jo Jo 

and (|79p . In 2D case, there is one extra degree of freedom which can be controlled in experiments, which is the 
aspect ratio of the observation region. If the aspect ratio of the observation region is very large, then the distribution 
function is essentially the same as in the one-dimensional case. Below we concentrate on the case with unity aspect 
ratio. Similarly to one-dimensional case, the integral 

/ / dx 1 dy 1 \\og({x 1 -x 2 ) 2 + { yi -y 2 ) 2 ) (A59) 
Jo Jo * 

can be evaluated analytically, which somewhat simplifies the numerical evaluation. However, the dimensions of 
integrals grow fast with m, and here we present results only up to m = 3 : 

0.805087n -0.0740n + 0.3342n 2 0.015n + 0.045n 2 + 0.241n 3 n , 1 . , k . 
Z 2n = l + + — 2 + — 3 + 0(— A ). (A60) 

The error in numerical coefficient due to integration is of the order of the last reported digit. Compared to ID case, 
convergence is faster, which is consistent with our general ideology that the role of fluctuations is larger in systems 
with lower dimensions. 



APPENDIX B: JACK POLYNOMIALS 

(p) 

The microcanonical partition functions Z 2r l c an be com puted using so-called Jack polynomials. The Jack polynomi- 
als belong to a class of symmetric polynomials [1281 . Il7l| and are a one-parametric generalization of Schur symmetric 
functions and generically can be defined as 

Jx(x;g) = ^2 v *A0)' ln n( x )i ( B1 ) 
where the monomial symmetric function 771^(2;) 

cr 

includes sum over all permutations a. Here are numerical coefficients and i>a.a = 1- The ordering in sum (|B1[) 
means the ordering on the set of partitions A = (Ai, A2, . . . , A n ), Ai > A2 > . . . > A„ > and /1 = (j/,±,fi2, . . . , //„), 
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A*i > A*2 > ■ • ■ > A*n > of natural numbers A and /i. The ordering /i < A here means that 53j=i Mi ^ f° r an y 

fe < n — 1. The partitions can be represented using Young diagrams as follows: we put Ai boxes in the first row, A2 
boxes in the second, and so on. The definition of Jack polynomials thus includes a sum over all Young diagrams for 
which the number of rows is smaller or equal to the number of variables n. The Jack polynomials have t he pr operty of 
orthogonality which allows to bring the microcanonical partition function into the following form (see [120l | for more 
details on derivation) 

7 (p) = 2 v^tt + 2g("-* + l)]r[p + Aj + ^(n-i + l)] 
2n nZ.il r[A, + l + ^(n-i)]r[p + A i + ^(n-i)] • [ > 

Here c n =T[n + 1]/T n [l/2K] and sum goes over Young diagrams labeled by integers Ai > A2 > . . . A n > 0. 

In particular, the Jack polynomials expressions give the following results for the lowest microcanonical partition 
functions 



Z ( P) = BjnC&jrq-i/g) 

sin7r(^ +p)T(l - & + p)T(l - & - p) 



and 



4 " T^) [ ^(1-^)^(1 + ^) >**" L > K' K>' iL+ 2K' L+ 2K>' L) 



DC 



}j + — +n,l + — +n},{2 + n,2 + n},l)], (B5) 

n=0 ^ ' 

where p F q is a generalized hypergeometric function jl72| . Apparently, for larger n, expressions for become more 
and more cumbersome and are difficult to use for real computations. 

Jack polynomials appear particularly in connection to Calogero-Sutherland model, i.e. the model which describes n 
particles (the same n as in formulas above) on a line i ntera cting via inverse square interaction. In the first quantized 
form the Hamiltonian of this model can be defined as jl28l | 
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j — 1 J l<2<j<n v J ' 

Wave functions for the excited states of this model are Jack polynomials. 



APPENDIX C: THERMODYNAMIC BETHE ANSATZ OF THE QUANTUM IMPURITY MODEL 

The mo del described by the partition function can be solved exactly by thermodynamic Bethe ansatz (TBA) 
173, 174], This solution is non perturbative in g. One should be careful about the correspondence between the 
perturbative expansion and the results from the non perturbative TBA calculations. The correct correspondence 
between Z(K,g) pert and the non perturbative result of the TBA is given by 

Z (K,g) pert = — t== gxp[Ftba {K, g)] , (CI) 
V 2iv 

where Ftba{K, g) is the free energy. It can be expressed in terms of the energies e(9) of elementary excitations which 
depend on rapidity 9 : 

/oo jq 2K 1 

.„ ^ ^[Ff^j Ml + exp(e +W) ). (C2) 



Here, 

and the variable x is the fugacity in the perturbative expansion 



L \2K> L \2K-l) 
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To clarify the issue further we briefly describe the formalism of TB A. T he TBA equations are parametrized by the 
Dynkin diagrams corresponding to the classical algebras A n , D n , E n [175| . Our case under consideration is related to 
the diagrams D n , which describes the sine-Gordon model. To each node of the diagram we associate the particle with 
the mass given by the Coxeter number of the corresponding diagram (which is equal to 2K — 2) shown in Fig. (fTT)) . 
These are denoted by /i + = /i_ (corresponding to the kink and antikink) and fij, j = 1...2K — 2, (note that we consider 
K to be of the form n/2, n = 3, 4, 5, ...) corresponding to breathers, bound states of kink and antikink. Explicitly, 

fij = 2/i+ sin(^^ — ) for breathers j = 1...2K - 2. (C4) 




FIG. 17: Dynkin diagram for the algebra D n . In the context of sine-Gordon problem nodes 1...2K — 2 correspond to breathers, 
whereas nodes ± correspond to kink/antikink. 

Due to a great simplification of the structure of scattering matrices in integrable theory jl75l | the free energy of the 
boundary sine-Gordon problem can be expressed entirely in terms of the energies e+ (see Eq. (|C2j) above). This can 
further be reduced as follows. Consider the TBA equations for the energies: 

/°° d6' 2K — 1 

.„ 2^ cosh[(2if-l)(^)] ^ ^ l0g(1 + (C5) 

where N rs is the incidence matrix of the corresponding Dynkin diagram. To clarify these notations consider few 
examples. 

For K = 3/2 we have kink, antikink and one breather. Therefore the diagram consists of the three nodes, say 
+, 1, — . The vertex + (kink) is connected to the node 1 and the vertex 1 is connected to the vertex — . Vertices + 
and — are disconnected. The incidence matrix is therefore 

(C6) 

(the number of rows is equal to the number of vertices, the number of columns is equal to the number of bonds 
between these vertices). Using this matrix we obtain the following set of coupled TBA equations: 

/oo inf o 
^ c o^-n '° gll+eM ' im> - <C7) 

/oo inf n 
.^ cosh^^O] 106 ^^^^ (C8) 

where we use the fact that e+ = e_. Now, from the second equation we can see that the £i(8)/2 gives precisely the 
function F from the Eq. (|C2|) . Therefore, for the case of K = 3/2, 




3 11 

Zo{ 2 , 9)pert = exp[-ei (a)], (C9) 



where t\ must be determined from the set of coupled equations (|C7|) . 

To continue, consider the case of K — 4/2. Here we have 2 breathers, kink, and antikink. The incidence matrix is 

(CIO) 
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The set of TBA equations is therefore 



' (">=/ ^ —nJt flA1 log(l + exp( e2 (fl'))), (Cll) 



2?r cosh[3(6» - i 


9')] 


tW 3 




2tt cosh[3(6» - i 


9')] 1 


dff 3 




3 2tt cosh[3(6> - 




dS 7 3 


l 



<*(*)=/ — ^„uro/ fl flA1 log(l + exp( e2 (gO)), (C12) 



+2 i- 0o ^ coah^-y)] l0g(1 + W))' ( C13 ) 

from which we conclude that the £2/3 gives the function F from the Eq. (|C2[) . Therefore, for the case of X = 4/2, 

4 11 

Zo{^,g) P ert = -^=exp[-e 2 (a)], (C14) 

where e 2 must be determined from the set of coupled equations (|C11|) . 
For K = 5/2 we can find that 

5 11 

Zo(^,g)pert = ^=exp[-(e 3 (a) -e x (a))], (C15) 

where £1,3 must be determined from the set of corresponding coupled equations (easy to write down explicitly). 
For K = 6/2 we have 

Zo{\,9) P ert = ^= ex P [i(e 4 - e 2 + F( Cl ))], (C16) 

and therefore it includes the function F(e\). This is not convenient from the point of view of numerics since it contains 
an additional integration. 

Continuing the above computations one can observe that for the values of K = 2k/2, k = 2,3, 4, ... we will always 
have an additional function F(ei) in the expression for the perturbative function Zo(g, K) pert . On the other hand 
for the values of K = (2k + l)/2, k = 1,2,3, ... one can always express the result for partition function in terms of 
energies e%. It seems that this is more convenient for numerics. The result for this case is the following: 

Zo(g,^^) P ert = -^=i=exp[i(^(-l) s e (2n _ 1) _ 2s )], (C17) 

~ v s— 

where the set of energies must be determined from the solution of coupled integral equations. 

It is useful to consider the asymptotic behavior. For a large 9 —> ±00 values of the integration parameter, the 
energies of breathers may be taken in the following form (see Eq. (|C4[l ): 

e s (a) = — ^- = 2sin[ lexpfal, as a — > 00. (C18) 

/Z+ — 2 

The general strategy is clear now: one can solve numerically the set of coupled integral equations for some values 
of K = (2n + l)/2 and compute the partition function using Eq. (|CT7l) . At this point the difficulty arises: how 
to unambiguously define the analytic continuation of Z p (K,g) into imaginary-^ domain. To avoid this problem it is 
better to use the alternative approach described in section ITVC 11 
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